基于剪切效應纖維梁單元的結構非線性有限元數(shù)值模擬*

李嘉鈺，陳夢成，王開心

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

引 言

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

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

1 考慮剪切效應的纖維模型

1.1 基本假設

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

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

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

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

3） 剪應力與剪應變均勻分布.

1.2 剛度矩陣

1.2.1 截面剛度矩陣

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

與其材料性質相關的本構關系矩陣為

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

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

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

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

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

1.2.2 彈性剛度矩陣

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

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

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

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

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

其中

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

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

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

1.2.3 幾何剛度矩陣

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

非線性有限元方程為

其中，BNL為非線性應變矩陣，ΔP為等效荷載.

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

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

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

式中G和的顯式為[4]

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

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

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

2.1 等效應力與等效應變

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

等效應力與等效應變存在以下關系：

其中，φ為已知的材料本構關系.

2.2 基于增量理論應力的更新

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

1） 根據已知的應力分量 σ和應變分量ε 計算出等效應力和 等效應變，計算應變分量增量 dε.

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

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

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

然后，計算應力分量增量 dσ：

最后，按第5）步更新應力.

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

其中

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

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

然后，計算應力分量增量 dσ：

最后，按第5）步更新應力.

5） 得到纖維單元的應力增量后，更新應力：

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

圖2 增量理論應力更新流程圖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，考慮結構自重的影響，將重力等效為柱頂?shù)募辛，同時，在頂部作用一水平集中荷載V，考慮幾何非線性和材料非線性，對此結構做水平力-側移的非線性全過程分析.

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

算例分析采用本文所給出的三維纖維梁單元建模計算，計算時將原結構劃分為10 個單元，單元的每個節(jié)點有6 個自由度；劃分纖維單元時，盡可能保證原鋼筋面積與纖維網格的面積相等，優(yōu)先考慮此條件；將截面劃分為31×31的正方形纖維網格，網格按照結構的實際情況分為鋼筋纖維和混凝土纖維，如圖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 混凝土本構模型

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

圖5 混凝土本構模型Fig.5 The constitutive model of concrete

其中

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

3.1.2 鋼筋本構模型

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

圖6 鋼筋雙折線本構模型Fig.6 The constitutive model of steel reinforcement

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

3.1.3 數(shù)值模擬結果

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

圖7 荷載-位移關系曲線Fig.7 The load-displacement curves

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

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

3.2 算例二

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

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

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

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

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

4 結 論

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

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

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