李嘉鈺，陳夢成，王開心

（華東交通大學 土木建筑學院,南昌 330013）

引 言

基于梁單元的有限元分析模型是近年來人們關(guān)注的熱點，使用不同梁單元理論的有限元模型對所分析結(jié)構(gòu)的計算精度和計算效率有顯著影響.早期的Euler-Bernoulli 經(jīng)典梁理論[1]應(yīng)用廣泛，在一般情況下可以得到比較滿意的結(jié)果，但由于未考慮剪切變形給截面帶來的影響，在梁高度和跨度相差不大的情況下，會有較大的誤差.為了考慮剪切帶來的影響，Timoshenko[2]于1921年提出了考慮剪切效應(yīng)的修正梁理論，該理論仍保留了平截面假定，但實際上梁截面的剪應(yīng)力與剪應(yīng)變呈拋物線分布，是不均勻的，這與平截面假定相矛盾，于是引入了截面修正系數(shù)，將截面上不均勻的剪應(yīng)力與剪應(yīng)變等效為均勻分布，這樣既能較好地考慮剪切的影響同時又簡化了計算模型.Mari 和Scordelis[3]提出了基于有限元剛度法的纖維梁模型，該模型基于Euler-Bernoulli 梁理論，以單元積分點截面處的力學性能代表整個單元的力學性能，使用線性插值構(gòu)造積分點的軸向位移場，用Hermite 插值函數(shù)構(gòu)造橫向位移場.胡鄭州等[4]在UL 列式下，根據(jù)連續(xù)性介質(zhì)力學和虛位移原理，并結(jié)合Timoshenko 梁理論，推導出了考慮剪切效應(yīng)的三維纖維梁單元模型.顧觀東[5]在Euler-Bernoulli 梁單元模型的基礎(chǔ)上，修正了塑性階段的平截面假定，并引入剪切位移，得到了更精確的梁單元模型.鄭欣怡[6]將Timoshenko 梁理論引入單軸本構(gòu)的彎剪耦合纖維模型中，這種改進后的纖維梁單元模型考慮了彎剪耦合效應(yīng)和剪切效應(yīng)對軸向變形的間接影響.Thai 等[7]利用纖維梁單元，對單向軸力和彎矩作用下的橋墩進行了雙非線性有限元分析.張傳超等[8]基于OPENSEES 軟件，對纖維梁單元在鋼管混凝土柱應(yīng)用方面的建模方法進行了研究.馬銀[9]基于纖維梁單元，對型鋼混凝土結(jié)構(gòu)進行建模分析，并進一步推導了考慮黏結(jié)滑移關(guān)系的型鋼混凝土梁柱單元模型.藺鵬臻等[10]基于ABAQUS 軟件，采用精細化的纖維梁單元模型，開發(fā)了鋼筋和混凝土纖維梁單元材料用戶子程序.由此可見，纖維梁單元在理論研究方面和工程應(yīng)用方面都得到了廣泛的應(yīng)用[11-13].

目前，基于Euler-Bernoulli 梁理論的纖維梁單元其纖維只能考慮單軸受力狀態(tài)，無法考慮剪切效應(yīng)；基于考慮剪切效應(yīng)的Timoshenko 梁理論，多以使用材料多維本構(gòu)模型來考慮材料非線性，需計算的參數(shù)較多，且運算復雜.為此，本文基于考慮剪切效應(yīng)的三維纖維梁單元有限元模型，推導了該纖維梁單元的截面剛度矩陣、單元彈性剛度矩陣和單元幾何剛度矩陣，同時補充了文獻[4]在推導插值函數(shù)矩陣時，矩陣元素N2,6的遺漏，并根據(jù)UL 列式法和彈塑性增量理論，考慮結(jié)構(gòu)的幾何非線性和材料非線性，旨在建立壓彎剪復雜應(yīng)力狀態(tài)下結(jié)構(gòu)非線性有限元分析理論，能在合理描述壓彎剪作用的同時得到準確的計算結(jié)果.最后，使用MATLAB 編制了相關(guān)計算程序，結(jié)合鋼筋混凝土和矩形鋼管混凝土的典型壓彎剪構(gòu)件進行有限元數(shù)值模擬，驗證了本文所建立理論的有效性、正確性和通用性.

1 考慮剪切效應(yīng)的纖維模型

1.1 基本假設(shè)

本文基于經(jīng)典纖維模型，并結(jié)合Timoshenko 梁理論，建立了考慮剪切效應(yīng)的纖維模型.纖維梁單元模型如圖1所示，并做如下假設(shè)：

圖1 纖維梁單元Fig.1 The fiber beam element

1） 滿足平截面假定，但截面不再與梁中軸線垂直，考慮剪切效應(yīng)導致的截面翹曲；

2） 鋼筋與混凝土充分黏結(jié)，無相對滑移，變形協(xié)調(diào)；

3） 剪應(yīng)力與剪應(yīng)變均勻分布.

1.2 剛度矩陣

1.2.1 截面剛度矩陣

對于截面內(nèi)任一纖維單元，其空間幾何變換矩陣為

與其材料性質(zhì)相關(guān)的本構(gòu)關(guān)系矩陣為

其中，Ei和Gi分別為第i個纖維單元的彈性模量和剪切模量.

通過積分運算，得到任一纖維單元的剛度：

其中，Ai為纖維單元的截面面積.

對截面內(nèi)的全部纖維單元進行積分，得到截面剛度矩陣：

將式(1)、(2)代入式(4)中，展開并進行積分運算，可以得到截面剛度矩陣的顯式為

1.2.2 彈性剛度矩陣

在得到截面剛度矩陣的基礎(chǔ)上，可采用位移形函數(shù)將結(jié)構(gòu)的位移與截面的變形聯(lián)系起來，從而推導三維纖維梁單元的剛度矩陣.每個梁單元共有2 個節(jié)點，每個節(jié)點有6 個自由度，單元節(jié)點位移向量為

其中，(u1,u2),(v1,v2),(w1,w2)分別為節(jié)點在x軸、y軸和z軸方向的線位移，(θ1x,θ2x)為節(jié)點的扭轉(zhuǎn)角位移，(θ1y,θ2y)和(θ1z,θ2z)分別為節(jié)點在xOz和xOy平面內(nèi)的轉(zhuǎn)動角位移.

單元軸線上任一點x處截面的位移用單元兩端節(jié)點位移可表示為

式中，N(x)為單元形函數(shù)矩陣.軸向位移及扭轉(zhuǎn)角位移均采用相同的線性位移模式，橫向彎曲位移均采用三次式的位移模式[14]，此時單元位移插值函數(shù)可寫為

單元形函數(shù)矩陣為

其中

對單元形函數(shù)矩陣進一步推導，可得到線性應(yīng)變矩陣BL，再由虛功原理和變分原理得單元彈性剛度矩陣為

做積分運算后，得到單元彈性剛度矩陣的顯式為

其中，單元彈性剛度矩陣為對稱矩陣，即

1.2.3 幾何剛度矩陣

在兩節(jié)點12 個自由度單元線彈性分析力學模型的基礎(chǔ)上，本文將建立單元的幾何非線性分析力學模型，并推導出單元的幾何剛度矩陣.

非線性有限元方程為

其中，BNL為非線性應(yīng)變矩陣，ΔP為等效荷載.

在式(12)中，左邊第一項積分為彈性剛度矩陣，通過第二項積分，可推導幾何剛度矩陣.

其中，ATσ=Mθ=MGΔue，A為非線性應(yīng)變與位移梯度矢量之間的轉(zhuǎn)換矩陣，M為應(yīng)力分量矩陣，θ為位移梯度矢量，為位移梯度矢量與單元節(jié)點位移矢量之間的轉(zhuǎn)換矩陣，故有G

因此，單元幾何剛度矩陣為

式中G和的顯式為[4]

其中，σ11為沿x軸的軸向應(yīng)力，τxy為垂直x軸沿y軸的剪應(yīng)力，τxz為垂直x軸沿z軸的剪應(yīng)力，τyz為垂直y軸沿z軸的剪應(yīng)力.

通過式(15)的積分運算可得單元幾何剛度矩陣，最終，單元剛度矩陣為彈性剛度矩陣和幾何剛度矩陣之和，即

2 復雜應(yīng)力狀態(tài)下纖維單元應(yīng)力的更新

2.1 等效應(yīng)力與等效應(yīng)變

本文將采用Mises 屈服準則判斷復雜應(yīng)力狀態(tài)下纖維單元的屈服狀態(tài)，并據(jù)此理論求等效應(yīng)力與等效應(yīng)變.式(17)和(18)為等效應(yīng)力與等效應(yīng)變的分量表示方法：

等效應(yīng)力與等效應(yīng)變存在以下關(guān)系：

其中，φ為已知的材料本構(gòu)關(guān)系.

2.2 基于增量理論應(yīng)力的更新

由于塑性本構(gòu)關(guān)系中的應(yīng)力和應(yīng)變不存在一一對應(yīng)的關(guān)系，一般只能建立應(yīng)力與應(yīng)變增量之間的關(guān)系，以增量形式來表示.依據(jù)此理論對纖維單元應(yīng)力狀態(tài)的更新步驟如下：

1） 根據(jù)已知的應(yīng)力分量 σ和應(yīng)變分量ε 計算出等效應(yīng)力和 等效應(yīng)變，計算應(yīng)變分量增量 dε.

2） 在Mises 屈服準則下，通過比較等效應(yīng)力與纖維單元的屈服應(yīng)力判斷是否屈服，若未屈服，按第3）步彈性理論更新應(yīng)力；若屈服，按第4）步塑性理論更新應(yīng)力.

3） 根據(jù)彈性理論，計算彈性Jacobi 矩陣De：

其中，E為纖維單元的彈性模量，μ為Poisson 比.

然后，計算應(yīng)力分量增量 dσ：

最后，按第5）步更新應(yīng)力.

4） 根據(jù)塑性理論，計算塑性Jacobi 矩陣Dep：

其中

式(23)中，G為纖維單元的剪切模量；為通過式(17)所求得的等效應(yīng)力；Ep為塑性模量，它與彈性模量E和切線模量Et之間的關(guān)系為Ep=EEt/(E?Et)，切線模量可通過對本構(gòu)關(guān)系的求導得到.

平均應(yīng)力、應(yīng)力偏量按式(24)計算：

然后，計算應(yīng)力分量增量 dσ：

最后，按第5）步更新應(yīng)力.

5） 得到纖維單元的應(yīng)力增量后，更新應(yīng)力：

更新材料應(yīng)力狀態(tài)的流程圖如圖2所示.

圖2 增量理論應(yīng)力更新流程圖Fig.2 The flow chart for stress updating with the incremental theory

3 算例分析

3.1 算例一

一正方形截面鋼筋混凝土柱[15]，如圖3所示，柱高1.65 m，截面尺寸為550 mm×550 mm，縱向受力鋼筋直徑為20 mm，等距分布，保護層厚度為40 mm，考慮結(jié)構(gòu)自重的影響，將重力等效為柱頂?shù)募辛，同時，在頂部作用一水平集中荷載V，考慮幾何非線性和材料非線性，對此結(jié)構(gòu)做水平力-側(cè)移的非線性全過程分析.

圖3 鋼筋混凝土柱（單位：mm）Fig.3 The reinforced concrete column (unit:mm)

算例分析采用本文所給出的三維纖維梁單元建模計算，計算時將原結(jié)構(gòu)劃分為10 個單元，單元的每個節(jié)點有6 個自由度；劃分纖維單元時，盡可能保證原鋼筋面積與纖維網(wǎng)格的面積相等，優(yōu)先考慮此條件；將截面劃分為31×31的正方形纖維網(wǎng)格，網(wǎng)格按照結(jié)構(gòu)的實際情況分為鋼筋纖維和混凝土纖維，如圖4所示，其余參數(shù)見表1.

表1 計算參數(shù)Table 1 Calculation parameters

圖4 纖維單元截面Fig.4 The fiber element section

本文迭代求解的方法采用杜修力等提出的位移控制新方法[16]，控制點為水平力施加點，位移增量控制步步長=1 mm，總位移增量步數(shù)n=25，迭代收斂準則采用位移收斂準則，允許誤差δ=10?3.

3.1.1 混凝土本構(gòu)模型

考慮到混凝土受到箍筋約束作用，故在數(shù)值模擬中，受拉區(qū)混凝土采用經(jīng)Scott 等[17]修正后的Kent-Park[18]單軸混凝土本構(gòu)模型，如圖5所示.需要說明的是，對于受拉區(qū)混凝土，經(jīng)試算表明，可以忽略其受拉貢獻，即取受拉混凝土纖維的拉應(yīng)力為0，這一簡化在不影響計算精度的前提下提高了計算效率.混凝土抗壓強度及相應(yīng)的混凝土峰值應(yīng)變和極限應(yīng)變均取文獻[15]中所給出的數(shù)值.混凝土受壓本構(gòu)模型具體表達式如下：

圖5 混凝土本構(gòu)模型Fig.5 The constitutive model of concrete

其中

式中，K為考慮箍筋約束效應(yīng)所引起的混凝土強度增大系數(shù)，根據(jù)本算例參數(shù)可計算得K=1.21；ρsv為箍筋的體積配筋率；fyh為箍筋的屈服強度；εu為混凝土極限壓應(yīng)變，本文取 εu=0.024 8；Zm為約束混凝土應(yīng)力應(yīng)變曲線下降段的斜率，經(jīng)式(29)計算，本文取Zm=13.66.

3.1.2 鋼筋本構(gòu)模型

鋼筋的本構(gòu)模型采用雙折線彈塑性本構(gòu)模型，如圖6所示，此本構(gòu)模型具體表達式如下：

圖6 鋼筋雙折線本構(gòu)模型Fig.6 The constitutive model of steel reinforcement

式中，Es為鋼筋的彈性模量；fy為鋼筋的屈服強度；εy為鋼筋的屈服應(yīng)變；k為鋼筋硬化段斜率，本文與文獻[15]相同，取k=0.01Es.

3.1.3 數(shù)值模擬結(jié)果

利用彈塑性增量理論結(jié)合考慮剪切效應(yīng)下的三維纖維梁單元模型，編制相應(yīng)的MATLAB 程序計算，將數(shù)值模擬結(jié)果與文獻[15]的數(shù)據(jù)進行對比分析，如圖7和表2所示.

圖7 荷載-位移關(guān)系曲線Fig.7 The load-displacement curves

表2 數(shù)值模擬結(jié)果Table 2 Numerical simulation results

由圖7和表2可見：本文的計算結(jié)果與文獻[15]結(jié)果吻合較好，水平荷載的全過程有限元數(shù)值模擬結(jié)果略低于文獻[15]的結(jié)果，其相對誤差均在8%以內(nèi)，結(jié)構(gòu)的彈性階段、彈塑性階段以及下降段均較為明顯.上述結(jié)果說明，本文所建立的復雜應(yīng)力狀態(tài)下考慮剪切效應(yīng)的纖維梁單元理論是正確的，依據(jù)該理論所編制的非線性有限元分析程序，能有效地解決實際問題.

3.2 算例二

一矩形鋼管混凝土壓彎剪構(gòu)件[19]，試件長度L=1 200 mm，截面高D=300 mm，截面寬B=200 mm，外鋼管厚度t=5.6 mm，含鋼率 α=0.1，所用混凝土強度等級為C60，鋼管為Q345 鋼，試件端部所受軸壓N=0.4Nu（Nu為該構(gòu)件的軸壓極限承載力），跨中強軸方向受水平荷載V.在數(shù)值模擬有限元分析模型中，試件兩端均為鉸接，鋼材的本構(gòu)模型采用雙折線模型，混凝土的本構(gòu)模型采用塑性損傷模型[20].計算簡圖如圖8所示，采用本文所給出的考慮剪切效應(yīng)的纖維梁單元模型建模，使用位移控制新方法[16]進行有限元非線性方程的迭代求解，計算水平荷載V與支座處轉(zhuǎn)角位移R(R=2Δ/L，Δ為外荷載V產(chǎn)生的水平位移）的全過程關(guān)系曲線.

圖8 矩形鋼管混凝土壓彎剪試件計算模型Fig.8 The calculation model for the rectangular CFST column under the compression-bending-shear loading condition

圖9為本文數(shù)值模擬結(jié)果與文獻[19]的V-R曲線.由圖可知，兩者在彈性階段吻合程度良好，在彈塑性階段和下降段，由于纖維梁單元模型不能模擬鋼管和混凝土的相互作用，本文的計算結(jié)果稍滯后于文獻[19]的結(jié)果，但偏差相對較小，從總體上來說，兩者結(jié)果較為吻合.此外，本文計算出的構(gòu)件所受最大橫向力為479.8 kN，文獻[19]的最大橫向力約為480.1 kN，說明本文的纖維梁單元模型能較好地模擬構(gòu)件的橫向極限承載力.

圖9 矩形鋼管混凝土試件V-R 曲線Fig.9 V-R curves of the rectangular CFST column

綜上所述，本文所建立的復雜應(yīng)力狀態(tài)下結(jié)構(gòu)非線性分析理論是正確的，證明了依據(jù)此理論所編制的非線性程序是可行、有效的.

4 結(jié) 論

1） 算例分析的結(jié)果表明，本文基于考慮剪切效應(yīng)的纖維梁單元，所建立的壓彎剪復雜應(yīng)力狀態(tài)下結(jié)構(gòu)非線性有限元分析理論是通用、可行和正確的.

2） 利用本文理論進行結(jié)構(gòu)的非線性全過程分析時，可較好地模擬出結(jié)構(gòu)的彈性階段、彈塑性階段和下降段，且特征點較為明顯.

3） 依據(jù)本文理論所編制的非線性有限元分析程序是有效、可靠的，可以較好地解決理論研究和工程應(yīng)用中所涉及的結(jié)構(gòu)非線性問題.