基于非線性屈曲材料模型的張弦桁架連續倒塌分析

鄭義凡，曾 濱，周 臻，許 慶

(1. 東南大學混凝土及預應力混凝土結構教育部重點實驗室，南京 210096；2. 中冶建筑研究總院有限公司，北京100088)

隨著社會經濟水平的提高和工程技術水平的高發展，大跨空間結構被廣泛運用于會展中心、體育館等大空間的公共建筑。在大跨空間結構應用廣泛的同時，結構連續倒塌事故也時有發生[1-3]，造成了慘重的人員傷亡和經濟損失，大跨結構連續倒塌問題不容忽視。

在目前的大跨空間結構有限元分析中，古莉等[4]為考慮桿件的屈曲效應，將每根桿件劃分為多個單元并施加初始幾何缺陷。增量動力分析法需要進行大量的有限元分析，采用此方法將增加模型的計算時間，分析效率低，難以適用于大規模的增量動力分析，也缺乏對張弦桁架抗連續倒塌性能的易損性評估。

在有限元動力分析研究中，張弦桁架在構件失效后因動力效應，結構產生震蕩，桿件單元在進入塑性階段后可能會經歷加載、卸載的循環作用，因此在材料模型中需考慮材料模型的強化問題和單元的滯回能力。隨動強化模型為常用的鋼材本構模型[5-6]，其假定加載過程中的后繼屈服面在應力空間作剛體移動而沒有轉動，初始屈服面大小、形狀和方向保持不變。試驗表明：后繼屈服面既有曲面大小的變化，也有中心位置的移動。為了彌補隨動強化模型的不足，DING 等[7]經過研究后提出了適用于梁柱單元的混合強化本構模型。

對于傳統張弦桁架結構連續倒塌分析，俞鋒等[8]、武嘯龍[9]和張超等[10]多集中于研究結構連續倒塌過程中某些桿件的軸力時程曲線以及某些節點的位移時程曲線，從而得出桿件軸力以及節點位移的變化規律，但缺乏對張弦桁架抗連續倒塌性能的評估，從而無法評價張弦桁架在不同工況下結構的抗連續倒塌性能。

為此，本文在梁柱單元修正混合強化本構模型的基礎上，結合梁柱單元非線性屈曲材料，提出了張弦桁架連續倒塌快速分析方法，在建立有限元模型時只需將桿件劃分為一個單元，可同時兼顧分析效率與精度。利用ANSYS/LS-DYNA 編制了提出材料模型的子程序，開展了基于IDA 的張弦桁架抗連續倒塌性能評估。最后，通過概率分析得出了張弦桁架在拉索失效工況下連續倒塌易損性曲線。

張弦桁架在構件失效后因動力效應，結構產生震蕩，桿件單元在進入塑性階段后可能會經歷加載、卸載的循環作用，因此在材料模型中需考慮材料模型的強化問題和單元的滯回能力。

大多數網殼結構的桿件為圓鋼管，而網殼結構空間效應較強，導致圓鋼管任意一點的應力僅包括軸向應力和扭轉產生的剪應力[11]。DING 等[7]根據以上結論提出了梁柱單元混合強化本構模型，其后繼屈服函數如下：

1 梁柱單元非線性屈曲材料模型

然而，對于單榀張弦桁架結構，其空間效應較弱，在復雜受力工況下桿件中存在不可忽略的彎矩。梁柱單元普通混合強化本構模型只存在扭轉產生的剪應力τxy和軸向應力 σx，因此需要引入第二個方向上的剪應力以修正該本構模型，即τxz，從而使單元應力狀態符合單榀張弦桁架結構以及其他梁柱結構的實際情況。

根據Zeigler 修正運動硬化法則，加載曲面沿聯結其中心和現時應力點的向量方向移動，此時后繼屈服函數可表示為：

式中：M為混合強化參數，取值范圍為[0, 1]，當M=1 時為等向強化，M=0 時為隨動強化，M取其它值時表示混合強化，一般取值為0～0.2，本文取0.2[12]； dλ為一非負的比例系數。

為保證桿件在劃分為一個單元時能夠考慮壓桿的屈曲效應，需開展材料模型的二次開發工作。圖1 為HILL 等[13]提出的非彈性后屈曲本構模型，HILL 模型的非線性屈曲段采用指數函數進行描述，但是函數中的兩個參數X1、X2被定為常數[14]，與桿件長細比無關，缺乏一定的合理性。其中，曲線卸載段終點為點3，坐標為(0.5, 0.5)[15]。

圖1 非彈性后屈曲本構模型Fig. 1 Inelastic post-buckling constitutive model

謝道清等[16]采用ANSYS/LS-DYNA 建立了鋼管桿件的有限元模型，按照一階特征值屈曲模態引入桿件初彎曲，最大撓度取v0=L/1000。通過對46 根桿件進行加載卸載，得到一系列樣本桿件拉壓滯回P-Δ曲線，進行歸一化滯回曲線分析，最后以三次多項式擬合曲線控制參數得到等效彈塑性滯回模型，如圖2。模型中系數 α1、 α2、 α3、α4、 γ1為關于長細比 λ的參數。

圖2 等效彈塑性滯回模型Fig. 2 Equivalent elastic-plastic hysteretic model

綜上所述，本文采用Hill 模型的受壓屈曲卸載段(即圖1 中點2～點3 段)、等效彈塑性滯回模型的受壓屈曲段(即圖2 中點3～點6 段)，基于梁柱單元修正混合強化本構模型，得到梁柱單元非線性屈曲材料模型，如圖3，并進行材料模型的二次開發。利用該材料模型子程序，可以實現張弦桁架連續倒塌的快速分析。

圖3 梁柱單元非線性屈曲材料模型Fig. 3 Beam-column element nonlinear buckling material model

2 梁單元算例分析

為了驗證所編寫的梁柱單元非線性屈曲材料模型子程序的正確性，對標準桿件算例進行了非線性屈曲材料模型桿件、初始幾何缺陷模型桿件的軸壓分析。施加幾何初始缺陷的桿件初彎曲最大撓度v0=L/1000，采用BEAM161 單元模擬，材料模型采用*MAT_PLASTIC_KINEMATIC。對于使用本文所提出的材料模型桿件，采用BEAM161單元模擬，沿桿長方向劃分1 個單元。加載按位移控制，速度為 10-3m/s。

桿件規格分別為 Φ89×4 mm(L=3 m， λ=100)、Φ140×4.5 mm(L=3 m， λ=62)，E=2.06×1011Pa，σy=235 MPa。算例計算結果如圖4 所示。

圖4 桿件 P-Δ曲線Fig. 4 P-Δ Curve of member

由圖4 算例結果可知，梁柱單元非線性屈曲材料模型和初始幾何缺陷模型計算結果有著較小的差距，其原因是梁柱單元非線性屈曲材料模型是基于梁柱單元修正混合強化本構模型編寫，而初始幾何缺陷模型基于隨動強化模型。對于λ=100的桿件，梁柱單元非線性屈曲材料模型峰值承載力為175.34 kN，初始幾何缺陷模型峰值承載力為174.21 kN，二者相差0.65%。Δ =0.02 m時，梁柱單元非線性屈曲材料模型承載力為48.75 kN，初始幾何缺陷模型承載力為52.50 kN，二者相差7.14%。對于 λ=62的桿件，梁柱單元非線性屈曲材料模型峰值承載力為410.39 kN，初始幾何缺陷模型峰值承載力為407.88 kN，二者相差0.61%。Δ=0.02 m時，梁柱單元非線性屈曲材料模型承載力為140.72 kN，初始幾何缺陷模型承載力為148.95 kN，二者相差5.53%。以 λ=100的桿件為例進行滯回分析，滯回模擬結果如圖5 所示。

圖5 桿件(λ=100)滯回曲線Fig. 5 Hysteresis curve of member (λ=100)

由以上兩個算例可知，本文所編寫的梁柱單元非線性屈曲材料模型子程序能夠較好地模擬受壓桿件的非線性屈曲效應，與施加初始幾何缺陷桿件的P-Δ曲線基本吻合。同時桿件單元擁有良好的滯回能力，其受壓屈曲卸載段的終點為(0.5Fy,0.5 Δy)，符合理論材料模型。

3 張弦桁架連續倒塌試驗有限元模型驗證

采用梁柱單元非線性屈曲材料模型對已有張弦桁架連續倒塌動力試驗[9]開展數值模擬，并與試驗結果對比分析。結構跨度6 m，矢高0.4 m，垂度0.4 m。構件截面規格如下：上弦桿 Φ 20×2 mm，下弦桿Φ 20×2 mm，腹桿 Φ8 mm，撐桿Φ 32×2.5 mm，拉索 Φ8 mm。采用有限元軟件ANSYS/LS-DYNA建立分析模型，桁架桿件單元采用BEAM161，拉索單元采用LINK167，采用初始應變法施加預應力。鋼材屈服應力 σy=345 MPa，彈性模量E=2.06×1011Pa，極限應變取0.2；拉索預應力T=2 kN，彈性模量E=1.90×1011Pa；荷載取P=1041 N，均勻施加在桁架上弦節點；初始失效構件為跨中拉索，構件失效時間tp=0.2 s。

本文取兩個測點進行試驗數據與數值模擬數據對比。試驗中，B1 測點測量結構滑動鉸支座的水平位移，B2 測點測量結構跨中下弦節點的豎向位移，測點布置如圖6 所示。t=0 s 時開啟斷索裝置時拉索斷裂，t=0.243 s 時B1 測點水平位移達到最大值，因此時程曲線的時間范圍選為0 s～0.243 s。試驗與數值模擬時程曲線對比如圖7 所示。

圖6 測點布置圖Fig. 6 Layout of Measuring points

圖7 對比結果表明：數值模擬與試驗位移曲線總體基本一致，桁架倒塌破壞形式類似，整體符合較好，說明基于本文材料子程序的張弦桁架連續倒塌快速分析方法具有可行性。

圖7 測點時程曲線對比Fig. 7 Comparison of measuring point time history curve

4 張弦桁架連續倒塌IDA 分析

目前針對張弦桁架的連續倒塌研究多集中于研究結構連續倒塌過程中某些桿件的軸力時程曲線以及某些節點的位移時程曲線，缺乏對張弦桁架抗連續倒塌性能的評估。增量動力分析法雖然多用于結構抗震性能評估，但可以借鑒其思想，提出適用于張弦桁架連續倒塌的增量動力分析法，并利用IDA 曲線進行張弦桁架抗連續倒塌性能評估。其主要分析思路如下：

1)選取跨中下弦節點撓跨比作為DM(Damage measure)，結構外荷載值作為IM(Intensity measure)。

2) 令IM 由0 逐步增加至結構極限荷載，進行備用荷載路徑法分析，得到一系列的DM 值，以此繪出單條DM-IM 曲線。本文中，選取拉索為失效構件。

3)由于單榀張弦桁架結構冗余度較低，抗連續倒塌能力較弱，因此定義構件失效時，結構即超越IO(Immediate occupancy)性態。參考JGJ7-2010《空間網格結構技術規程》[17]，定義撓跨比達到1/250 時，結構達到CP (Collapse prevention)性態；定義曲 線 α 值 ( α為 荷載放大系數)達到αmax的80%時，結構達到GI (Global instability)性態。

4.1 張弦桁架結構信息

本文以單榀張弦桁架為研究對象，張弦桁架跨度80.75 m，具體結構尺寸及半結構部分桿件編號見圖8。其中，B1 為跨中下弦節點，A1～A10 為上弦桿，A11～A18 為下弦桿。上部桁架結構采用倒三角空間形式，每一節間由四角錐基本單元構成。柱頂設置固定鉸支座與上部拱結構鉸接，柱底剛接，于結構跨中上弦桿兩段設置水平側向約束。構件截面規格如下：上弦桿 Φ165×10 mm，下弦桿 Φ215×16 mm，腹桿 Φ108×8 mm，撐桿Φ180×8 mm，柱帽桿 Φ325×18 mm，柱Φ700×30 mm，拉索2- 127Φ7 mm。采用有限元軟件ANSYS/LS-DYNA 建立分析模型，桁架桿件單元采用BEAM161，拉索單元采用LINK167，采用初始應變法施加預應力。桿件單元材料采用梁柱單元非線性屈曲材料模型，索單元材料采用*MAT_CABLE_DISCRETE_BEAM。鋼材屈服應力σy=235 MPa，彈性模量E=2.06×1011Pa，極限應變取0.2。拉索預應力T=400 kN，彈性模量E=1.90×1011Pa。基準荷載取P=5 kN，定義為 α荷載放大系數，外荷載 α×P均勻施加在桁架上弦節點。

圖8 張弦桁架立面示意圖 /mFig. 8 Elevation diagram of cable-stayed arch-truss

4.2 張弦桁架連續倒塌快速分析方法計算耗時

傳統的張弦桁架連續倒塌分析方法在建立有限元模型時通常將桿件劃分為多個單元，而本文所提出的快速分析方法只需將桿件劃分為一個單元。現分別建立三種有限元模型，M-1 模型采用快速分析方法，M-2、M-3 模型采用傳統分析方法，分別將桿件劃分為4 單元和6 單元。有限元計算耗時如表1 所示。

表1 有限元計算耗時表Table 1 Finite element calculation time

由表1 知，快速分析方法相比于傳統分析方法，其計算耗時相比于4 單元縮短了39.89%，相比于6 單元縮短了57.36%。因此，在開展增量動力分析，涉及計算模型較多時，快速分析方法將極大地減少計算耗時。同時，在進行有限元建模時，無需將桿件劃分為多個單元并施加初始幾何缺陷，提升了建模效率。

4.3 拉索失效工況下連續倒塌分析

α=1.0、1.2、1.4、1.6、1.8、1.9 時的B1 節點豎向位移時程曲線如圖9 所示， α=1.9 時上弦桿中A1 桿～A5 桿和下弦桿中A11 桿～A15 桿軸力時程曲線如圖10 所示。

圖9 跨中節點豎向位移時程曲線Fig. 9 Vertical displacement time history curve of mid-span node

圖10 桿件軸力時程曲線Fig. 10 Axial force time history curve of members

由圖9 知，當 α<1.9 時，B1 節點豎向位移隨荷載增加緩慢；當 α=1.9 時，跨中節點豎向位移曲線不再收斂，結構發生連續倒塌。

根據梁柱單元非線性屈曲材料模型，上弦桿屈曲荷載為-1039 kN，下弦桿屈服荷載為2350 kN。t=0 s 時，拉索失效。由圖10 知，上弦桿、下弦桿軸力迅速增大。當 α=1.9、t=0.28 s 時，A5 桿軸力達到-1039 kN，桿件發生屈曲，軸力下降，而A1 桿～A4 桿 仍 處 于 彈 性 階 段；A11 桿～A15 桿軸力均小于2350 kN，未屈服。由于A5 桿屈曲，導致A1 桿～A4 桿軸力下降，進而導致A11 桿～A15桿軸力下降。由以上分析知，對于本算例，當拉索失效時，位于桁架跨中的上弦桿首先受壓屈曲，繼而導致周圍桿件軸力下降，結構整體承載力降低，發生連續倒塌。為提高結構的抗連續倒塌能力，跨中上弦桿需要加強。跨中附近的下弦桿軸力接近屈曲荷載，因此，同樣需要加強以防止跨中下弦桿屈服。

4.4 桿件截面加強的張弦桁架IDA 分析

由4.3 節分析知，在拉索失效工況下， α =1.9、t=0.28 s 時，桁架跨中附近的上弦桿受壓屈曲，跨中附近的下弦桿處于彈性階段。因此，加強的桿件如表2 所示，S-1 為原結構，S-4 加強了所有屈曲上弦桿。上弦桿截面加強至 Φ205×10 mm，下弦桿截面加強至 Φ255×16 mm。

表2 加強桿件表Table 2 Reinforced members

以B1 節點撓跨比( λ ×10-2, λ=Δy/L×100 ,Δy為跨中下弦節點豎向位移)作為結構性能參數，四種結構的IDA 曲線如圖11 所示。四種結構在拉索失效工況下均經歷了類似的破壞階段，在曲線起始段結構存在明顯的彈性范圍；隨著 α的增大，曲線達到拐點，B1 節點豎向位移迅速增大，曲線斜率減小，并伴隨小幅度波動，結構進入彈塑性階段；當 α達到最大值時，結構發生連續倒塌。

圖11 IDA 曲線Fig. 11 IDA Curve

四種結構達到不同性態所對應的 α值如表3 所示。由圖11 知，對于 αmax，S-2 無提升，S-3、S-4提升16.67%；對于CP 性態，S-2 無提升，S-3 提升33.3%，S-4 提升83.3%；對于GI 性態，S-2 無提升，S-3、S-4 提升16.67%。由以上分析知，S-2相較于S-1 的抗連續倒塌能力提升較小， S-3、S-4相較于S-1 的抗連續倒塌能力提升較大，因此，加強部分重要桿件能夠提升張弦桁架在拉索失效工況下的抗連續倒塌能力。

表3 結構性態表Table 3 States of structures

以S-1 為例探究A5 桿軸力在外荷載逐步增大過程中的變化規律，S-1 結構A5 桿分別在α=0.4、0.6、0.8 時的軸力時程圖如圖12 所示。當α=0.4 時，A5 桿軸力在結構處于拉索失效后的振動狀態中，其最大軸力為-928 kN，桿件未屈曲；當 α=0.6 時，A5 桿最大軸力接近上弦桿屈曲軸力，此時圖11 中S-1 結構的IDA 曲線達到拐點，結構整體進入彈塑性階段；當 α=0.8 時，A5 桿軸力時程曲線的波動形勢不再呈圓滑的正弦狀，且結構穩定后A5 桿軸力小于 α=0.6 的A5 桿軸力，結合圖3 的梁柱單元非線性屈曲材料模型，表明在拉索失效后A5 桿件經歷了屈曲后卸載、加載的循環作用，導致桿件剛度下降，繼而導致結構整體剛度下降，IDA 曲線斜率減小。

圖12 S-1 結構A5 桿軸力時程曲線Fig. 12 Axial force time history curve of A5 member, S-1

參考結構地震易損性分析，進行張弦桁架連續倒塌易損性分析。假設 α 和 λ之間服從指數分布，其對數關系為：

通過IDA 分析數據的回歸分析，可以得到式(13)中的A和B。對于S-1，A=1.3704，B=0.3985；對于S-2，A=1.4076，B=0.4608；對于S-3，A=1.4753，B=0.2913；對于S-4，A=1.8072，B=-0.2056。

張弦桁架連續倒塌易損性計算模型為：

5 張弦桁架連續倒塌易損性分析

5.1 概率倒塌需求模型

5.2 拉索失效工況下連續倒塌易損性曲線

將式(14)代入式(13)并結合4.4 節中的性態點判斷即可得到張弦桁架在拉索失效工況下的連續倒塌易損性曲線，如圖13 所示。由圖13 知，在拉索失效工況下，三種結構超越CP 性態的概率最大，超越GI 性態稍有減小，即隨著張弦桁架結構破壞程度的增大，各性態水準的超越概率減小。除S-2 結構外，對于CP 性態和GI 性態，隨著加強桿件數量的增多，結構的超越概率逐漸減小。

圖13 拉索失效工況下連續倒塌易損性曲線Fig. 13 Progressive collapse fragility curve under cable failure conditions

S-2 結構相較于S-1 結構，CP 性態超越概率相差較小。且由圖11 知，當 α=80%αmax時，S-2結構的 λ小于S-1 結構的 λ，因此導致S-2 結構的GI 性態超越概率反而大于S-1 結構的GI 性態超越概率。表明當加強較少數量的重要桿件數量時，無法達到提高張弦桁架結構抗連續倒塌能力的目標。

對比不同大小荷載放大系數 α下，四種結構超越CP 性態和GI 性態的概率大小，如表4 所示。在 α=1.0、2.0、3.0 時，S-3 結構相比于S-1 結構對GI 性態的超越概率分別減小57.14%、14.07%、2.26%，S-4 結構相比于S-1 結構對GI 性態的超越概率分別減小98.84%、59.67%、14.04%。表明隨著荷載的增大，加強重要桿件對張弦桁架抗連續倒塌能力的提高程度逐漸下降；在相同荷載的條件下，加強更多數量的重要桿件能夠給張弦桁架抗連續倒塌能力帶來更高程度的提升。

表4 四種結構超越概率對比Table 4 Comparison of fragility probability between four structures

6 結論

本文提出了一種張弦桁架連續倒塌的快速分析方法，并針對試驗和工程算例開展了分析，結論如下：

(1) 基于梁柱單元修正混合強化本構模型，結合非彈性后屈曲本構模型和等效彈塑性滯回模型，編寫了ANSYS/LS-DYNA 梁柱單元非線性屈曲材料模型子程序，并通過與試驗結果的對比，驗證了快速分析方法的可行性。

(2) 針對張弦桁架結構算例，采用本文方法進行了單榀張弦桁架在拉索失效工況下的連續倒塌分析。統計了三種有限元模型的計算耗時，采用本文方法的有限元模型計算耗時至少縮短39.89%，極大減少了計算耗時。

(3) 通過原結構與三個加強結構的IDA 曲線對比分析，結果表明：在拉索失效工況下，加強跨中上弦桿、跨中下弦桿能夠提高張弦桁架的抗連續倒塌能力，極限承載力最多提高16.67%，達到CP 性態、GI 性態所對應的荷載放大系數也均有所提高；當荷載達到一定值時，拉索失效后跨中上弦桿經歷了屈曲后卸載、加載的循環過程，導致桿件剛度下降，繼而導致結構整體剛度下降，IDA 曲線出現拐點。

(4) 開展了張弦桁架連續倒塌的易損性分析，結果表明：基于快速分析方法可有效評估張弦桁架的綜合抗連續倒塌能力。