基于非线性屈曲材料模型的张弦桁架连续倒塌分析

郑义凡，曾 滨，周 臻，许 庆

(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) 开展了张弦桁架连续倒塌的易损性分析，结果表明：基于快速分析方法可有效评估张弦桁架的综合抗连续倒塌能力。