一種組合式屈曲單元及其在輸電塔結構倒塌仿真分析中的應用

付 興，鐘璽峰，李宏男,2，朱 宇

(1. 大連理工大學海岸和近海工程國家重點實驗室，遼寧，大連 116024；2. 沈陽建筑大學土木工程學院，遼寧，沈陽 110168)

鋼結構由于其構造簡單和性能優(yōu)良的特點，在輸電塔、通信塔、工業(yè)廠房等結構中得到了廣泛應用。這類鋼結構構件長細比較大，在強風、地震等極端荷載作用下構件的破壞模式以屈曲為主。特別是輸電塔結構，屬于典型的空間鋼結構，在強風作用下受壓構件很容易發(fā)生屈曲破壞[1-5]，因此，如何準確評估這類結構的屈曲行為及承載能力是結構抗災設計的關鍵。

我國地處季風性氣候區(qū)，尤其是東南沿海地區(qū)，臺風和特大暴雨等天氣時常出現(xiàn)。在風雨荷載共同作用下，高聳柔性的大跨越輸電塔線體系經常發(fā)生振動疲勞損傷和連續(xù)性倒塌破壞[6-7]。針對輸電塔的連續(xù)倒塌破壞，國內外學者開展了大量的研究工作。Albermani 等[8]用非線性方法分析并模擬了輸電塔的倒塌過程，進而開展了輸電塔倒塌的足尺實驗。Patil 等[9]在輸電塔倒塌模擬中考慮了材料和幾何非線性的影響，與線彈性分析結果對比表明，考慮各種非線性因素后桿件內力有較大提升?；诮Y構多尺度分析技術，王鳳陽[10]考慮幾何非線性和材料非線性的影響，對輸電塔倒塌過程中變形較大的桿件建立了精細殼單元模型，模擬了下?lián)舯┝飨螺旊娝牡顾^程以及桿件內力的變化。姚旦等[11]采用向量式有限元方法分析了某輸電塔結構的倒塌過程。

國內外學者提出了很多力學模型用于計算和模擬軸心受壓構件的屈曲行為和滯回特性，可分為三大類：彈塑性有限元法、塑性鉸法和經驗模型法[12]。彈塑性有限元法是用材料的本構關系來模擬構件的滯回性能。Jin 和El-Tawil[13]考慮了沿截面高度方向的塑性發(fā)展，將有限元法用于鋼框架的動力非線性分析；Salawdeh 和Goggins[14]考慮構件疲勞斷裂，提出一種新的支撐單元模型。彈塑性有限元法的優(yōu)點在于精度高、通用性強，缺點是占用計算機內存較大，顯式算法增量步長不能太大，使得計算效率較低，且隱式算法迭代需要設置合理的參數(shù)。塑性鉸法是基于支撐的受力變形特點，用彈性梁單元和塑性鉸組成簡化的支撐模型。Dicleli 和Calik[15]推導出支撐滯回曲線的力與位移的關系，考慮增長效應和鮑辛格效應，并將數(shù)值模擬結果與實驗對比，驗證了該方法的準確性；Takeuchi 和Matsui[16]研究了支撐在隨機循環(huán)載荷下的斷裂狀態(tài)。塑性鉸法優(yōu)點在于模型簡單、計算效率高，缺點是精度不如有限元法。經驗模型法基于實驗數(shù)據(jù)建立簡化的桿件軸向力-位移模型來模擬支撐的滯回行為。Haroun 和Shepherd[17]采用經驗模型法對鋼框架進行動力時程分析，研究了支撐的非線性行為對結構抗震性能的影響；劉慶志等[18]提出一種新的鋼支撐經驗模型法，提高了其模擬精度和計算效率。經驗模型法優(yōu)點是計算速度快、收斂性好，缺點是模型所需的參數(shù)不易確定。

在輸電塔結構倒塌仿真分析中，大部分學者采用彈塑性有限元法模擬構件的力學特性。為了更高效準確地模擬輸電塔構件的非線性屈曲行為，結合塑性鉸法的優(yōu)點，本文首先在第一節(jié)提出了一種組合式屈曲單元，并對其構造和工作原理進行了詳細的闡述，第二節(jié)介紹了自編程序采用的靜動力方程求解方法，第三節(jié)對MATLAB 編寫的仿真程序進行了全面驗證，最后第四節(jié)梳理全文工作。

1 組合式屈曲單元構造原理

細長構件在兩端鉸接約束下，塑性鉸一般出現(xiàn)在構件中點；在兩端剛接約束下，塑性鉸多產生于中點和兩端；一端剛接一端鉸接的約束情況下，塑性鉸通常出現(xiàn)在固定端和距離固定端0.7 倍構件長度的位置[19]。

為高效模擬輸電塔構件屈曲行為的線性及非線性力學特性，本文提出了一種組合式屈曲單元。該單元在兩端剛接情況下塑性鉸位置會出現(xiàn)在中間及兩端，所以該組合式單元由兩個傳統(tǒng)彈性梁單元組成，并在兩個傳統(tǒng)彈性梁單元端部布置轉動彈簧，在連接兩個彈性梁單元的中間位置布置拉壓彈簧和轉動彈簧，拉壓彈簧單元及轉動彈簧單元均采用理想彈塑性力-位移模型。同時為了更加有效地模擬真實構件的屈曲行為，設置參數(shù)偏心距，考慮初始缺陷帶來的影響。下面以兩端剛接的細長構件為例，詳細介紹該組合式屈曲單元構造，其他約束形式的組合式屈曲單元原理相同、構造方法一致，本文就不再詳細介紹。

兩端剛接的組合式屈曲單元自由度編號如圖1所示。從單元中間塑性鉸位置截斷，拆分為兩個彈性梁單元，左邊梁單元的平動自由度為1～3、16～18，其轉動自由度為13～15、19～21；右邊梁單元的平動自由度為7～9、25～27，其轉動自由度為22～24、28～30。1～3、4～6 分別為節(jié)點n1的平動自由度和轉動自由度，7～9、10～12 分別為節(jié)點n6的平動自由度和轉動自由度?？紤]構件的初始缺陷，梁單元中間節(jié)點偏心距e取1/m單元長度。

圖1 組合式屈曲單元構造及自由度編號Fig. 1 Structure and serial numbers of degrees of freedom for combined buckling element

利用三個拉壓彈簧模擬滑動鉸，分別連接圖1中的16～18 與25～27 編號的平動自由度，當連接的自由度發(fā)生相對位移時，拉壓彈簧將產生反力限制其位移。當反力達到拉壓彈簧的屈服強度Fs時，拉壓彈簧屈服，屈服方向為原梁單元上、下兩節(jié)點連線方向，屈服后該方向的剛度為零，其余兩個方向的剛度不變。因為構件的長細比較大，不考慮橫向剪切破壞，軸向屈服強度采用如下公式計算：

式中：fy為鋼材屈服強度；A為截面面積。拉壓彈簧彈性剛度矩陣為：

設彈性梁單元軸向剛度為k，彈簧剛度ke等于彈性梁單元軸向剛度的n倍，即：

式中：E為單元彈性模量；l為組合式屈曲單元總長度。引入拉壓彈簧后的單元整體軸向剛度誤差為：

軸向剛度的折減比例為1/(n+1)，為保證軸向剛度，誤差在一定范圍內，n應取盡量大的數(shù)。由于剛度在位移計算中做除數(shù)，考慮計算機計算精度限制問題，除數(shù)不可取過大；同時，在動力分析中，若n取過大的數(shù)將會提高結構的最大頻率，從而需要更小的計算步長，會影響計算效率。基于上述考慮，本文n取1000，可保證軸向剛度誤差在0.1%以內。

將九個轉動彈簧分為三組，模擬三個塑性鉸，分別連接圖1 中4～6 與13～15、10～12 與22～24、19～21 與28～30 轉角自由度，當連接的自由度發(fā)生相對位移時，轉動彈簧將產生反力限制其位移。當反力達到轉動彈簧的屈服彎矩Ms時，轉動彈簧屈服，屈服方向為相應自由度的位移矢量和之差，屈服后該方向的剛度為零，其余兩個方向不變。轉動彈簧的屈服彎矩Ms等于考慮軸力時的最大截面合力矩。以角鋼為例，構件全截面進入塑性后的受力示意圖如圖2 所示，其中散點陰影區(qū)域為受拉區(qū)，面積為A1，斜線陰影區(qū)域為受壓區(qū)，面積為A2。截面合力為fy(A1-A2)，等于此刻單元的軸力，由以上條件計算得到的截面合力矩即為轉動彈簧的屈服彎矩Ms。同拉壓彈簧，轉動彈簧彈性剛度也取構件彎曲剛度的1000 倍。

圖2 角鋼構件屈服截面示意圖Fig. 2 Schematic diagram of yield section for angle steel member

在進行連續(xù)倒塌仿真時，失效構件的變形通常非常大，因此，本文桿件受拉失效準則定義如下：

式中：δ 為構件伸長率；l0和l分別為組合式屈曲單元的初始長度和變形后長度；δT,cr為鋼材的臨界伸長率，取20%。當構件伸長率超過臨界值，認為構件失效，將刪除中間的拉壓彈簧和彎曲彈簧，此時組合式屈曲單元分成兩個獨立的彈性梁單元。

該組合式屈曲單元采用拉壓彈簧和轉動彈簧實現(xiàn)了塑性鉸功能，簡化了塑性鉸法的編程復雜度、降低了剛度矩陣維度、大幅縮短了計算時長。

2 靜動力分析方法介紹

2.1 靜力分析

靜力分析即分析結構在給定靜力載荷作用下的響應，并求得結構位移、應力和應變等數(shù)據(jù)。本文自編求解程序考慮結構的幾何非線性，采用更新拉格朗日列式法建立靜力平衡方程[20]：

式中：m為單元個數(shù)；T為單元從整體坐標系到局部坐標系的轉換矩陣； Δu為整體坐標系下的節(jié)點位移增量；F為外荷載列陣；為單元彈性剛度矩陣；為單元幾何剛度矩陣；為大位移剛度矩陣，是由大位移引起的結構剛度變化，采用增量求解時可忽略[20]。

為避免上述方程隱式迭代計算中由于壓桿屈曲失穩(wěn)帶來的剛度突變，進而導致程序不收斂的問題，本文求解靜力平衡方程采用顯式方法，計算時根據(jù)非線性程度將靜力平衡計算分為N步，第i步的荷載為Fi=F×i/N，并根據(jù)靜力平衡方程進行每一步的平衡計算。

2.2 動力分析

在采用本文提出的組合式屈曲單元進行動力分析時，結構整體質量矩陣、剛度矩陣及荷載向量的組裝方法均和傳統(tǒng)有限元理論一致，不再贅述。采用Rayleigh 公式生成阻尼矩陣時，為避免結構進入非線性后，剛度矩陣變化帶來的阻尼矩陣不穩(wěn)定，本文Rayleigh 阻尼只考慮結構整體質量矩陣。

求解動力平衡方程時采用隱式-顯式混合計算方法，隱式算法采用Newmark-β 法，顯式算法采用中心差分法。在時程分析前期，無材料非線性行為(未出現(xiàn)塑性鉸)，幾何非線性影響也較小(位移較小)，隱式算法所需迭代的次數(shù)較少，此時，采用隱式算法的計算效率明顯高于顯式算法。當結構出現(xiàn)強非線性或是桿件發(fā)生失穩(wěn)破壞導致結構整體倒塌時，隱式算法無法收斂，此時，應該采用顯式算法。因此，本文編制的求解器首先采用隱式算法進行動力時程分析，當隱式算法不收斂時，將隱式算法最后一個收斂步的結構運動狀態(tài)作為顯式算法初始狀態(tài)，程序內核調用顯式算法繼續(xù)計算。

3 組合式屈曲單元及求解程序驗證

本文提出的組合式屈曲單元需通過自編求解程序進行結構靜動力分析，采用MATLAB 程序開發(fā)了基于組合式屈曲單元的輸電塔結構靜動力仿真平臺。為說明自編程序的精度及新單元的適用性，本章將與通用有限元分析軟件ABAQUS 及實驗結果進行全面對比驗證。

3.1 數(shù)值驗證

為驗證自編程序在強幾何非線性、線彈性狀態(tài)下，自編程序的剛度矩陣、質量矩陣和有限元求解方法的準確性，本章節(jié)將與解析解、通用有限元分析軟件進行一系列對比。

3.1.1 解析解對比

用自編程序及ABAQUS 軟件分別建立材料、截面、長度均相同的懸臂梁模型，截面為圓環(huán)，兩個模型桿件均劃分10 個單元，且忽略自重影響，具體結構參數(shù)詳見表1。在梁端施加大小相同的遞增荷載，作用方向垂直于桿件，圖3 詳細對比了ABAQUS 軟件、自編程序及解析解的荷載-位移結果，可以看出，不同計算方法的坐標點完全重合，說明提出的組合式屈曲單元及自編求解程序具備很高的求解精度。

表1 懸臂梁結構信息Table 1 Information of cantilever beam structure

圖3 懸臂梁荷載-位移結果對比Fig. 3 Comparison of load-displacement of cantilever beam

3.1.2 幾何非線性對比

用自編程序及ABAQUS 軟件分別建立材料、截面、長度均相同的懸臂梁模型，截面為正方形，邊長0.05 m，均劃分20 個單元，并在梁端施加相同大小的集中力計算結構響應，具體結構及荷載參數(shù)詳見表2。將荷載分10 步施加在該懸臂梁端部，作用方向垂直于桿件，且與正方形邊長平行，因為，考慮幾何非線性后結構響應無解析解，所以，圖4僅對比了ABAQUS 軟件及自編程序的荷載-位移結果，可以看出，坐標點完全重合，在荷載最大值時ABAQUS 計算的懸臂梁末端位移為0.8035 m，自編程序對應結果為0.8029 m，相對誤差僅為0.07%，說明提出的組合式屈曲單元及自編程序對幾何非線性具備較高的仿真精度。

表2 懸臂梁結構及荷載信息Table 2 Information of cantilever beam structure and load

圖4 大變形下懸臂梁荷載-位移對比Fig. 4 Comparison of load-displacement of cantilever beam considering large deformation

3.1.3 動力分析對比

采用自編程序及ABAQUS 軟件分別建立材料、截面、尺寸相同的兩層框架，如圖5 所示，結構底層中心為坐標原點，各節(jié)點坐標見表3，同高度各節(jié)點x、y坐標絕對值一致，所有構件參數(shù)一致，具體見表4。結構前十階自振頻率對比詳見表5，最大相對誤差只有0.41%，與ABAQUS 軟件計算結果非常接近。

表3 框架結構節(jié)點坐標Table 3 Node coordinates of frame structure

表4 框架結構構件參數(shù)Table 4 Element parameters of frame structure

表5 框架結構前十階頻率對比Table 5 Comparison of first ten frequencies of frame structure

圖5 框架結構有限元模型Fig. 5 Finite element model of frame structure

在x軸方向對該框架結構進行地震動時程分析，忽略結構材料和幾何非線性行為，采用遷安波，峰值加速度調幅為1 m/s2。9 號節(jié)點x方向位移時程的對比如圖6 所示，兩條曲線幾乎重合，最大誤差僅為0.12%，說明自編程序的動力仿真求解器具有較高的計算精度。

圖6 框架結構位移時程曲線對比Fig. 6 Comparison of time-history curves of displacement for frame structure

3.2 實驗驗證

3.2.1 構件滯回特性對比

以Black 等[21]的實驗數(shù)據(jù)為基礎，選取4 根鋼支撐進行數(shù)值模擬，鋼支撐具體參數(shù)如表6 所示，Brace1 和Brace2 為兩端鉸接，Brace3 和Brace4為一端剛接、一端鉸接。

表6 鋼支撐參數(shù)[21]Table 6 Parameters of steel brace[21]

圖7 為自編程序和實驗結果的滯回曲線對比，紅色虛線為自編程序計算結果，黑色實線為實驗結果。從圖中可以看出，自編程序計算結果與實驗結果較為吻合，本文提出的組合式屈曲單元能有效地模擬鋼支撐受壓屈曲后的負剛度現(xiàn)象。

圖7 構件滯回曲線對比Fig. 7 Comparison of hysteresis curves of members

圖中曲線受拉段的高度不相等，這是由于模擬受拉屈服的拉壓彈簧采用理想彈塑性力-位移模型，而實際構件受拉屈服后進入強化階段，抗拉強度會有所提高。細長構件在外荷載作用下通常發(fā)生屈曲失穩(wěn)，因而，采用理想彈塑性力-位移模型可以很好地模擬該類結構的倒塌破壞行為。

偏心距作為組合式屈曲單元的重要參數(shù)，對計算結果有較大影響，現(xiàn)以表6 中Brace1 桿件為例，采用不同的偏心距參數(shù)計算最大受壓承載力，詳見表7。從表中可以看出，當偏心距為L/300 時，其最大受壓承載力與實驗結果最接近，同時偏心距越小其最大受壓承載力越大。

表7 偏心距對承載力的影響Table 7 Effect of eccentricity distance on strength

3.2.2 輸電塔足尺實驗對比

采用Fu 等[22]文獻中的足尺輸電塔模型，用自編程序建立輸電塔仿真模型，進行靜力推覆分析和動力倒塌分析，并與實驗結果對比。實驗塔加載時，定義設計荷載為參考荷載，然后，從零開始分級加載，每級加載參考荷載的5%，加載至參考荷載的115%時實驗塔倒塌，動力倒塌分析加載時，在15 s 內將荷載從零線性遞增至參考荷載的115%，在15 s～20 s 保持荷載不變，在20 s～30 s將荷載線性遞增至參考荷載的120%，倒塌過程如圖8 所示。

圖8 實驗塔倒塌過程[22]Fig. 8 Collapse process of experimental tower[22]

自編程序靜力推覆分析加載的參考荷載和實驗塔完全相同，每級加載參考荷載的1%，加載至參考荷載的116%時倒塌，倒塌破壞位置如圖9 所示，圖中圓點表示塑性鉸位置，塑性鉸均出現(xiàn)于塔腿處，與實驗塔破壞位置相同。自編程序模擬的承載能力與實驗結果非常接近，進一步驗證了本文提出的組合式屈曲單元的準確性和適用性。對比實驗和仿真的連續(xù)倒塌動畫和倒塌破壞位置可知，各時刻的倒塌狀態(tài)非常相近，說明了本文開發(fā)的動力仿真平臺在連續(xù)倒塌模擬方面具有很高的精度。

圖9 推覆分析中輸電塔塑性鉸位置Fig. 9 Position of plastic hinges of transmission tower in pushover analysis

為說明本文提出的組合式屈曲單元及自編程序在計算效率方面的優(yōu)勢，現(xiàn)以上述輸電塔結構為例，自編程序生成的剛度矩陣維度為24 519×24 519，通用有限元軟件為了保證計算精度，將每個桿件劃分為5 個單元，生成剛度矩陣維度為383 721×383 721。同時以100%的參考荷載為例，自編程序計算時長3.3254 s，通用有限元軟件計算時長25.501 25 s，電腦CPU 型號為i5-8500，內存16 GB。由此可見，本文提出的組合式屈曲單元及自編程序可大幅度提高計算效率，且具備較高計算精度。

4 結論

本文基于有限元理論，提出了一種組合式屈曲單元模擬細長桿件的屈曲失穩(wěn)現(xiàn)象，囊括了三種常見的約束類型：

(1)該單元采用拉壓彈簧模擬單元軸向塑性伸長，采用轉動彈簧模擬單元塑性彎曲，實現(xiàn)了塑性鉸功能，簡化了塑性鉸法的編程復雜度，降低了剛度矩陣維度，大幅縮短了計算時長。

(2)基于MATLAB 程序編寫了該單元的靜動力仿真程序，可以開展工程結構的靜力推覆分析和動態(tài)倒塌仿真分析。

(3)通過和解析解、ABAQUS 商用軟件、構件實驗、輸電塔足尺實驗對比，全面驗證了該組合式屈曲單元的準確性及自編仿真程序的可靠性。