基于Hillinger-Reissner變分的混合單元的開發(fā)與應用

丁 陽，賴程鋼，李 寧，李忠獻

基于Hillinger-Reissner變分的混合單元的開發(fā)與應用

丁 陽1,2，賴程鋼1，李 寧1,2，李忠獻1,2

(1. 天津大學建筑工程學院，天津 300072；
2. 天津大學濱海土木工程結構與安全教育部重點實驗室，天津 300072)

基于Hellinger-Reissner變分建立了一種用于幾何和材料非線性分析的混合單元，采用位移形函數(shù)和考慮二階效應的內力形函數(shù)對截面位移和截面內力進行插值，建立混合控制微分方程；通過靜力凝聚消除單元節(jié)點力未知量得到單元剛度矩陣和單元內力．單元狀態(tài)確定過程中，截面層次的平衡方程和單元層次的協(xié)調方程均通過引入非線性迭代算法以消除殘余誤差，從而減少結構層次的迭代次數(shù)．混合單元結合纖維截面模型用于鋼管混凝土(CFST)構件的數(shù)值分析，結果表明：相對于剛度法和柔度法單元，基于Hellinger-Reissner變分的混合單元可以更加準確地反映構件的幾何和材料非線性，非線性迭代算法用于單元狀態(tài)確定具有良好的計算效率和數(shù)值穩(wěn)定性．在此基礎上對影響CFST構件幾何非線性的主要參數(shù)進行了分析．

Hellinger-Reissner變分；混合單元；非線性迭代；幾何非線性；材料非線性；鋼管混凝土(CFST)

隨著工程結構復雜程度增大，結構非線性分析日漸發(fā)揮重要功能，廣泛應用于工程和研究領域．近年來，如何提高梁柱單元在結構幾何、材料非線性問題求解中的分析精度和計算效率一直是研究者關注的熱點問題之一．

基于纖維截面模型的分布塑性梁柱單元可準確地反映軸力和彎矩共同作用下材料的非線性，具有良好的分析精度和計算效率．然而，對于大變形狀態(tài)下的幾何非線性，一直未能得到很好的解決．

聶建國等[1]基于有限元程序MSC.MARC開發(fā)了用于組合結構地震反應分析的剛度法纖維梁單元；王文達等[2]基于非線性纖維梁-柱理論建立了分布塑性有限元模型，用于鋼管混凝土(concrete filled steel tube，CFST)框架結構全過程分析．傳統(tǒng)的剛度法單元采用三次橫向位移和線性軸向位移假定，適用于線性曲率和軸向常應變位移場，無法準確描述構件在強非線性狀態(tài)下的真實位移場分布形式．通過增加單元數(shù)量或者采用更高階形函數(shù)[3-4]來逼近真實位移場可以提高計算精度，但將影響單元計算效率．

相對于剛度法單元，柔度法單元以節(jié)點力作為基本未知量，沿構件全長插值得到截面內力，嚴格滿足單元平衡方程．但是現(xiàn)有的有限元程序中柔度法單元均為幾何線性單元，基于小變形假定建立單元控制方程，無法考慮構件在大變形下的幾何非線性．針對這一問題，陳滔等[5-6]結合二階柔度法和纖維模型梁柱單元，在內力形函數(shù)中引入橫向位移考慮構件二階效應．文獻[7-9]基于Hillinger-Reissner變分提出了考慮幾何非線性的分析方法，以單元節(jié)點位移和節(jié)點力作為基本未知量通過位移形函數(shù)和內力形函數(shù)插值得到截面位移和截面內力，最終通過靜力凝聚消除節(jié)點力未知量得到單元剛度矩陣．該方法基于兩場變分能夠較為準確地反映構件在大變形狀態(tài)下的受力狀態(tài)，并且相對柔度法單元在數(shù)值計算上具有更好的穩(wěn)定性．

本文在已有研究基礎上，基于Hellinger-Reissner變分建立混合單元有限元方程，將非線性迭代算法用于截面層次和單元層次消除殘余誤差，從而減少結構層次迭代次數(shù)．將混合單元結合纖維截面模型嵌入開源有限元程序OpenSees用于鋼管混凝土構件非線性分析，以充分驗證混合單元以及非線性迭代算法用于非線性分析的有效性．基于混合單元分析了構件長細比和荷載偏心比對鋼管混凝土構件幾何非線性的影響規(guī)律．

1 基于Hillinger-Reissner變分的混合單元

圖1給出了混合單元有限元方程建立的流程.混合單元以節(jié)點位移和節(jié)點力作為基本未知量，分別通過位移形函數(shù)和內力形函數(shù)插值得到截面位移和截面內力，由截面位移得到的截面應變與由截面內力得到的截面應變之間應滿足協(xié)調方程；截面內力和外荷載之間應滿足單元平衡方程；最終由單元平衡方程和協(xié)調方程基于Hellinger-Reissner變分建立單元剛度矩陣．單元狀態(tài)確定過程中，由節(jié)點力插值得到的截面內力與由截面應變積分得到的截面內力尚需滿足截面平衡方程．

混合單元基于以下基本假定：①不考慮剪切和扭轉變形的影響；②單元截面應變狀態(tài)符合平截面假定；③單元外荷載只考慮桿端集中荷載；④單元軸向變形滿足小變形假定．

1.1 單元力與變形

混合單元基本未知量包括單元節(jié)點位移向量U和節(jié)點力向量Q，即

式中：ΔL為桿端相對軸向變形；θzi、θzj為繞z軸的桿端轉角；θyi、θyj為繞y軸的桿端轉角；Nx為桿端軸力；Mzi、Mzj為繞z軸的桿端彎矩；Myi、Myj為繞y軸的桿端彎矩．

單元軸向位移和橫向位移分別采用線性形函數(shù)和三次Hermitian形函數(shù)，截面位移向量u可以表示為

式中：u(x)為截面軸向位移；v(x)、w(x)分別為截面2個主軸方向上的橫向位移；Nu為單元位移形函數(shù)矩陣．

為描述單元應變與位移的非線性關系，由單元截面位移向量u得到的截面應變向量d采用Green-Lagrange應變張量表示為

式中：ε為截面軸向應變；κz、κy為繞z、y軸的曲率．

單元內力采用常軸力和線性彎矩并考慮二階效應，單元截面內力向量q可表示為

式中：N(x)為截面軸力；Mz(x)、My(x)為繞z、y軸的截面彎矩；Nq為單元內力形函數(shù)矩陣．

式中L為單元長度．

1.2 截面平衡方程

基于纖維截面廣義截面內力與變形關系由截面變形向量dq積分得到截面內力向量qsum，與由節(jié)點力向量Q插值得到的截面內力向量q應滿足截面平衡方程

式中Dq為截面不平衡力．

1.3 單元剛度矩陣

基于虛位移原理，混合單元平衡方程的虛功方程可表示為

式中：δdTq為截面內力向量q所做的虛功；-δUTP為外荷載向量P所做的虛功．

基于纖維截面廣義截面內力與變形關系由截面內力向量q得到單元截面應變向量dq，與由單元截面位移向量u得到的截面應變向量d應滿足變形協(xié)調性，即

基于Hellinger-Reissner變分[10]，式(8)和式(9)可表示為

經(jīng)分部積分得到

式中E、C分別為單元求解過程中的不平衡力和殘余變形向量．

利用泰勒級數(shù)展開式可將非線性方程(12)～(13)線性化，具體過程在文獻[9]中已有介紹，這里不再詳述．

式中：ΔQ為單元節(jié)點力增量；ΔU為單元節(jié)點位移增量；ΔP為單元外力增量；N為單元軸力．

OpenSees是基于剛度法編制的有限元程序，為了實現(xiàn)將混合單元嵌入該程序，由式(14)、式(15)通過靜力凝聚將單元節(jié)點力未知量ΔQ消除得到

式中：Ke為混合單元剛度矩陣；Qe為單元內力向量．

1.4 單元狀態(tài)確定

混合單元需要對截面平衡方程、單元協(xié)調方程及單元平衡方程進行求解，單元狀態(tài)確定過程較為復雜，文獻[11]介紹了4種用于混合單元狀態(tài)確定的算法．Denavit等[9]采用線性算法通過線性求解將截面層次和單元層次的殘余誤差轉化為不平衡力在結構層次進行迭代求解．該算法需要在結構層次進行多次迭代以消除截面不平衡力和單元殘余變形的累積誤差，需要程序儲存大量單元矩陣．由于只在結構層次進行迭代求解，當單元層次殘余誤差過大時可能導致計算不收斂．

非線性迭代算法在結構層次每一次迭代過程中都會分別在截面層次和單元層次通過迭代將對應的殘余誤差消除，從而減少結構層次的迭代次數(shù)，程序只需儲存少量的單元矩陣．本文基于OpenSees采用C++語言編制程序將非線性迭代算法用于混合單元狀態(tài)確定，具體計算過程如圖2所示．

圖2 非線性迭代算法Fig.2 Nonlinear iteration algorithm

2 算例分析

2.1 材料本構關系

本文圓鋼管混凝土應力-應變曲線(見圖3)采用Denavit等[9]改進的Chang-Mander模型，混凝土應力-應變關系為

圖3 圓鋼管混凝土應力-應變曲線Fig.3 Stress-strain relationship curve of circular CFST

本文方鋼管混凝土應力-應變曲線(見圖4)在混凝土達到峰值應力前參考韓林海[12]模型考慮方鋼管對混凝土的約束效應，峰值應力之后參考文獻[13]考慮構件截面尺寸、鋼管和混凝土強度對混凝土應變軟化及殘余應力的影響．其應力-應變關系為

圖4 方鋼管混凝土應力-應變曲線Fig.4 Stress-strain relationship curve of square CFST

核心混凝土峰值應力后應變軟化曲線斜率Kc表示為

式中：B/t為方鋼管寬厚比；Es為鋼材彈性模量．

核心混凝土殘余應力frc和對應的壓應變εrc可表示為

2.2 軸壓試驗

Han[14]進行了一組細長圓鋼管混凝土柱軸壓試驗，本文選取其中2個試件進行分析，其試驗參數(shù)見表1，初始缺陷取L/5,000[15]．根據(jù)對稱性取一半結構采用1個單元模擬，單元設置4個積分點，截面離散為16個鋼纖維和32個混凝土纖維．混合單元的計算結果與試驗結果對比如圖5所示．可見，計算結果與試驗結果吻合良好，說明混合單元能較好地反映細長鋼管混凝土軸壓柱的極限承載力和軟化行為．

表1 軸壓試驗試件參數(shù)Tab.1 Test parameters of specimen under axial load

圖5 軸壓柱混合單元計算結果與試驗結果對比Fig.5 Comparison of mixed element computational results and test results of axial compression columns

2.3 偏壓試驗

Grauers[16]進行了一組方鋼管混凝土柱偏壓試驗，整個加載過程隨著荷載的增大試件在軸力和彎矩的共同作用下逐漸表現(xiàn)出材料非線性和幾何非線性.本文分別采用混合單元、剛度法單元和柔度法單元對其中7個方鋼管混凝土柱進行分析，試件S-2和S-10基本試驗參數(shù)見表2．

偏壓柱3種單元的計算結果與試驗結果對比如圖6所示．可見，整個加載過程中混合單元計算結果與試驗結果吻合良好，較為準確地反映了試件的初始剛度、極限承載力以及軟化剛度．基于剛度法單元、柔度法單元的計算結果在加載初期與試驗結果吻合較好，隨著試件變形增大計算結果與試驗結果誤差逐漸增大，均高估了試件極限承載能力．

表2 偏壓試驗試件參數(shù)Tab.2 Test parameters of specimen under eccentric load

圖6 偏壓柱3種單元的計算結果與試驗結果對比Fig.6Comparison of computational results and test results of three kinds of elements of eccentric compression columns

表3給出了7個偏壓試件的數(shù)值計算結果與試驗結果誤差．從表中可以看出，采用1個混合單元模擬即可得到較為精確的結果；剛度法單元和柔度法單元計算結果與試驗結果具有較大誤差，需要成倍增加單元數(shù)量才能減小計算結果與試驗結果誤差，嚴重降低了計算效率．

表3 偏壓構件數(shù)值計算結果與試驗結果誤差Tab.3 Errors of computational and test results of eccentric compression member

3 參數(shù)對構件幾何非線性的影響

以上分析說明不考慮幾何非線性會過高估計鋼管混凝土構件的承載力．本文針對試件S-10分析了參數(shù)對構件幾何非線性的影響，參數(shù)包括長細比λ和荷載偏心比e/B，其中λ＝40～140，e/B＝0.1～0.5．

圖7所示為長細比λ＝140、荷載偏心比e/B＝0.1的計算結果對比．可見，隨著變形的增大不考慮幾何非線性明顯高估構件承載能力，達到極限承載力60%左右時不考慮幾何非線性對構件承載力約高估10%，達到構件極限承載力時約高估16%．極限承載力之后相對誤差隨變形增大逐漸減小并趨于穩(wěn)定．

圖7 幾何非線性參數(shù)對構件承載力的影響Fig.7 Influence of geometric nonlinear parameters on bearing capacity of member

圖8 所示為構件長細比λ與荷載偏心比e/B對幾何非線性的影響．可見，荷載偏心比不變時，相對誤差與構件長細比大致呈線性增長關系，且小偏心比情況下相對誤差增長更快．構件長細比較小(λ≤80)時，極限承載力相對誤差隨荷載偏心比變化較為平緩，構件長細比較大(λ＞80)時，極限承載力相對誤差隨著荷載偏心比增大線性減小，說明大長細比情況下荷載偏心比對幾何非線性影響更明顯．

圖8 長細比與荷載偏心比對幾何非線性的影響Fig.8Influence of slenderness ratio and load eccentricity ratio on geometric nonlinearity

4 結 論

(1)建立了一種基于Hellinger-Reissner變分的混合單元，采用非線性迭代算法確定單元狀態(tài)，并將單元模型嵌入OpenSees用于非線性分析．

(2)采用非線性迭代算法在單元層次和截面層次將殘余變形消除，減少了結構層次的迭代次數(shù)，具有良好的穩(wěn)定性．相對于剛度法和柔度法單元，單個構件采用1個混合單元即可較為準確地反映構件強非線性狀態(tài)下的內力和變形，具有更好的分析精度和計算效率．

(3)不考慮幾何非線性對鋼管混凝土構件承載力大約高估3%～16%．長細比對幾何非線性有明顯的影響，長細比越大幾何非線性越明顯．構件長細比較小(λ≤80)時，荷載偏心比對幾何非線性影響較?。畼嫾L細比較大(λ＞80)時，隨著荷載偏心比增大，構件幾何非線性減弱．

［1］ 聶建國，陶慕軒. 采用纖維梁單元分析鋼-混凝土組合結構地震反應的原理[J]. 建筑結構學報，2011，32(10)：1-10.

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-10(in Chinese).

［2］ 王文達，韓林海. 鋼管混凝土框架結構力學性能非線性分析[J]. 建筑結構學報，2008，29(6)：75-86.

Wang Wenda，Han Linhai. Nonlinear finite elementanalysis on mechanical performance of concrete filled steel tubular frame structure[J]. Journal of Building Structures，2008，29(6)：75-86(in Chinese).

［3］ 夏擁軍，陸念力. 梁桿結構二階效應分析的一種新型梁單元[J]. 工程力學，2007，24(7)：39-43.

Xia Yongjun，Lu Nianli. A new beam element for second-order effect analysis of beam structures [J]. Engineering Mechanics，2007，24(7)：39-43(in Chinese).

［4］ 許紅勝，周緒紅，舒興平. 空間鋼框架幾何非線性分析的一種新單元[J]. 工程力學，2003，20(4)：39-44.

Xu Hongsheng，Zhou Xuhong，Shu Xingping. A new element for geometric nonlinear analysis of threedimensional steel frames[J]. Engineering Mechanics，2003，20(4)：39-44(in Chinese).

［5］ 陳 滔，黃宗明. 基于有限單元柔度法的材料與幾何雙重非線性空間梁柱單元[J]. 計算力學學報，2006，23(5)：524-528.

Chen Tao，Huang Zongming. Material and geometrically nonlinear spatial beam-column elements based on the finite element flexibility method[J]. Journal of Computational Mechanics，2006，23(5)：524-528(in Chinese).

［6］ 陳 滔，黃宗明. 基于有限單元柔度法和剛度法的幾何非線性空間梁柱單元比較研究[J]. 工程力學，2005，22(3)：31-38.

Chen Tao，Huang Zongming. Comparison between flexibility-based and stiffness-based geometrically nonlinear beam-column elements[J]. Engineering Mechanics，2005，22(3)：31-38(in Chinese).

［7］ Alemdar B N. Distributed Plasticity Analysis of Steel Building Structural Systems[D]. Atlanta：Department of Civil and Environmental Engineering，Georgia Institute of Technology，2001.

［8］ Tort C，Hajjar J F. Mixed finite element for threedimensional nonlinear dynamic analysis of rectangular concrete-filled steel tube beam-columns[J]. ASCE Journal of Structural Engineering，2010，136：1329-1339.

［9］ Denavit M D，Hajjar J F. Nonlinear seismic analysis of circular concrete-filled steel tube members and frames [J]. Journal of Structural Engineering，138(9)：1089-1098.

［10］ Zienkiewicz O C，Taylor R L. The Finite Element Method for Solid and Structural Mechanics[M]. 6th ed. Butterworth-Heinemann Ltd，Elsevier，2005.

［11］ Nukala P K，White D. Variationally consistent state determination algorithms for nonlinear mixed beam finite elements[J]. Computer Methods in Applied Mechanics Engineering，2004，193：3647-3666.

［12］ 韓林海. 鋼管混凝土結構：理論與實踐[M]. 2版. 北京：科學出版社，2007.

Han Linhai. Concrete Filled Steel Tubular Structures：Theory and Practice[M]. 2nd ed. Beijing：Science Press，2007(in Chinese).

［13］ Tort C，Hajjar J F. Mixed finite element approach for performance-based design of rectangular concrete-filled steel tube frames[C]//Composite Construction in Steel and Concrete Composite，Ⅵ-Proceedings of the 2008 Composite Construction in Steel and Concrete Conference. Tabemash，CO，USA，2011：591-603.

［14］ Han Linhai. Test on concrete filled steel tubular columns with high slenderness ratio[J]. Advances in Structural Engineering，2000，3(4)：337-344.

［15］ An Yufeng，Han Linhai，Zhao Xiaoling. Behavior and design calculations on very slender thin-walled CFST columns[J]. Thin-Walled Structures，2012，53：161-175.

［16］ Grauers M. Composite Columns with Hollow Steel Section Filled with High Strength Concrete [D]. Goteborg：Department of Structural Engineering and Mechanics，Chalmers University of Technology，1993.

(責任編輯：樊素英)

Development and Application of Mixed Element Based on Hellinger-Reissner Variation

Ding Yang1,2，Lai Chenggang1，Li Ning1,2，Li Zhongxian1,2
(1. School of Civil Engineering，Tianjin University，Tianjin 300072，China；2. Key Laboratory of Coast Civil Structure Safety of Ministry of Education，Tianjin University，Tianjin 300072，China)

This paper proposed a beam-column mixed element based on Hellinger-Reissner variation for geometric and material nonlinear analysis. A governing differential formulation of mixed element is established using displacement shape function and force shape function including second-order effects for interpolation of section deformation and innerforce，respectively．The element stiffness matrix and internal forces are derived by eliminating the force unknown through static condensation. During the element state determination，nonlinear iteration algorithm is adopted for residual error of section equilibrium and element compatibility，which reduces the numerical cost of structural equilibrium iteration. This mixed element based fiber section is used for numerical studies of concrete filled steel tube(CFST) members，and the results indicate good agreement. Relative to the displacement-based and flexibility-based elements，the mixed element based Hellinger-Reissner variation was more reasonable in predicting geometric and material nonlinearity of members，and the nonlinear iteration algorithm carried out for element state determination is of good efficiency and stability. Finally，a parametric study for geometric nonlinearity of CFST was performed based on the mixed element model.

Hellinger-Reissner variation；mixed element；nonlinear iteration；geometric nonlinearity；material nonlinearity；concrete filled steel tube(CFST)

O344.3

A

0493-2137(2015)05-0422-07

10.11784/tdxbz201310064

2013-10-25；

2013-11-20.

國家重點基礎研究發(fā)展計劃(973計劃)資助項目(2011CB013606，2011CB013603)；國家自然科學基金紫霞湖項目(51238007，

51378341，51178306)；天津市自然科學基金資助項目(13JCZDJC35200，13JCQNJC07200).

丁 陽（1966— ），女，博士，教授.

李 寧，neallee@tju.edu.cn.

時間：2013-12-31. 網(wǎng)絡出版地址：http://www.cnki.net/kcms/doi/10.11784/tdxbz201310064.html.