一种组合式屈曲单元及其在输电塔结构倒塌仿真分析中的应用

付 兴，钟玺峰，李宏男,2，朱 宇

(1. 大连理工大学海岸和近海工程国家重点实验室，辽宁，大连 116024；2. 沈阳建筑大学土木工程学院，辽宁，沈阳 110168)

钢结构由于其构造简单和性能优良的特点，在输电塔、通信塔、工业厂房等结构中得到了广泛应用。这类钢结构构件长细比较大，在强风、地震等极端荷载作用下构件的破坏模式以屈曲为主。特别是输电塔结构，属于典型的空间钢结构，在强风作用下受压构件很容易发生屈曲破坏[1-5]，因此，如何准确评估这类结构的屈曲行为及承载能力是结构抗灾设计的关键。

我国地处季风性气候区，尤其是东南沿海地区，台风和特大暴雨等天气时常出现。在风雨荷载共同作用下，高耸柔性的大跨越输电塔线体系经常发生振动疲劳损伤和连续性倒塌破坏[6-7]。针对输电塔的连续倒塌破坏，国内外学者开展了大量的研究工作。Albermani 等[8]用非线性方法分析并模拟了输电塔的倒塌过程，进而开展了输电塔倒塌的足尺实验。Patil 等[9]在输电塔倒塌模拟中考虑了材料和几何非线性的影响，与线弹性分析结果对比表明，考虑各种非线性因素后杆件内力有较大提升。基于结构多尺度分析技术，王凤阳[10]考虑几何非线性和材料非线性的影响，对输电塔倒塌过程中变形较大的杆件建立了精细壳单元模型，模拟了下击暴流下输电塔的倒塌过程以及杆件内力的变化。姚旦等[11]采用向量式有限元方法分析了某输电塔结构的倒塌过程。

国内外学者提出了很多力学模型用于计算和模拟轴心受压构件的屈曲行为和滞回特性，可分为三大类：弹塑性有限元法、塑性铰法和经验模型法[12]。弹塑性有限元法是用材料的本构关系来模拟构件的滞回性能。Jin 和El-Tawil[13]考虑了沿截面高度方向的塑性发展，将有限元法用于钢框架的动力非线性分析；Salawdeh 和Goggins[14]考虑构件疲劳断裂，提出一种新的支撑单元模型。弹塑性有限元法的优点在于精度高、通用性强，缺点是占用计算机内存较大，显式算法增量步长不能太大，使得计算效率较低，且隐式算法迭代需要设置合理的参数。塑性铰法是基于支撑的受力变形特点，用弹性梁单元和塑性铰组成简化的支撑模型。Dicleli 和Calik[15]推导出支撑滞回曲线的力与位移的关系，考虑增长效应和鲍辛格效应，并将数值模拟结果与实验对比，验证了该方法的准确性；Takeuchi 和Matsui[16]研究了支撑在随机循环载荷下的断裂状态。塑性铰法优点在于模型简单、计算效率高，缺点是精度不如有限元法。经验模型法基于实验数据建立简化的杆件轴向力-位移模型来模拟支撑的滞回行为。Haroun 和Shepherd[17]采用经验模型法对钢框架进行动力时程分析，研究了支撑的非线性行为对结构抗震性能的影响；刘庆志等[18]提出一种新的钢支撑经验模型法，提高了其模拟精度和计算效率。经验模型法优点是计算速度快、收敛性好，缺点是模型所需的参数不易确定。

在输电塔结构倒塌仿真分析中，大部分学者采用弹塑性有限元法模拟构件的力学特性。为了更高效准确地模拟输电塔构件的非线性屈曲行为，结合塑性铰法的优点，本文首先在第一节提出了一种组合式屈曲单元，并对其构造和工作原理进行了详细的阐述，第二节介绍了自编程序采用的静动力方程求解方法，第三节对MATLAB 编写的仿真程序进行了全面验证，最后第四节梳理全文工作。

1 组合式屈曲单元构造原理

细长构件在两端铰接约束下，塑性铰一般出现在构件中点；在两端刚接约束下，塑性铰多产生于中点和两端；一端刚接一端铰接的约束情况下，塑性铰通常出现在固定端和距离固定端0.7 倍构件长度的位置[19]。

为高效模拟输电塔构件屈曲行为的线性及非线性力学特性，本文提出了一种组合式屈曲单元。该单元在两端刚接情况下塑性铰位置会出现在中间及两端，所以该组合式单元由两个传统弹性梁单元组成，并在两个传统弹性梁单元端部布置转动弹簧，在连接两个弹性梁单元的中间位置布置拉压弹簧和转动弹簧，拉压弹簧单元及转动弹簧单元均采用理想弹塑性力-位移模型。同时为了更加有效地模拟真实构件的屈曲行为，设置参数偏心距，考虑初始缺陷带来的影响。下面以两端刚接的细长构件为例，详细介绍该组合式屈曲单元构造，其他约束形式的组合式屈曲单元原理相同、构造方法一致，本文就不再详细介绍。

两端刚接的组合式屈曲单元自由度编号如图1所示。从单元中间塑性铰位置截断，拆分为两个弹性梁单元，左边梁单元的平动自由度为1～3、16～18，其转动自由度为13～15、19～21；右边梁单元的平动自由度为7～9、25～27，其转动自由度为22～24、28～30。1～3、4～6 分别为节点n1的平动自由度和转动自由度，7～9、10～12 分别为节点n6的平动自由度和转动自由度。考虑构件的初始缺陷，梁单元中间节点偏心距e取1/m单元长度。

图1 组合式屈曲单元构造及自由度编号Fig. 1 Structure and serial numbers of degrees of freedom for combined buckling element

利用三个拉压弹簧模拟滑动铰，分别连接图1中的16～18 与25～27 编号的平动自由度，当连接的自由度发生相对位移时，拉压弹簧将产生反力限制其位移。当反力达到拉压弹簧的屈服强度Fs时，拉压弹簧屈服，屈服方向为原梁单元上、下两节点连线方向，屈服后该方向的刚度为零，其余两个方向的刚度不变。因为构件的长细比较大，不考虑横向剪切破坏，轴向屈服强度采用如下公式计算：

式中：fy为钢材屈服强度；A为截面面积。拉压弹簧弹性刚度矩阵为：

设弹性梁单元轴向刚度为k，弹簧刚度ke等于弹性梁单元轴向刚度的n倍，即：

式中：E为单元弹性模量；l为组合式屈曲单元总长度。引入拉压弹簧后的单元整体轴向刚度误差为：

轴向刚度的折减比例为1/(n+1)，为保证轴向刚度，误差在一定范围内，n应取尽量大的数。由于刚度在位移计算中做除数，考虑计算机计算精度限制问题，除数不可取过大；同时，在动力分析中，若n取过大的数将会提高结构的最大频率，从而需要更小的计算步长，会影响计算效率。基于上述考虑，本文n取1000，可保证轴向刚度误差在0.1%以内。

将九个转动弹簧分为三组，模拟三个塑性铰，分别连接图1 中4～6 与13～15、10～12 与22～24、19～21 与28～30 转角自由度，当连接的自由度发生相对位移时，转动弹簧将产生反力限制其位移。当反力达到转动弹簧的屈服弯矩Ms时，转动弹簧屈服，屈服方向为相应自由度的位移矢量和之差，屈服后该方向的刚度为零，其余两个方向不变。转动弹簧的屈服弯矩Ms等于考虑轴力时的最大截面合力矩。以角钢为例，构件全截面进入塑性后的受力示意图如图2 所示，其中散点阴影区域为受拉区，面积为A1，斜线阴影区域为受压区，面积为A2。截面合力为fy(A1-A2)，等于此刻单元的轴力，由以上条件计算得到的截面合力矩即为转动弹簧的屈服弯矩Ms。同拉压弹簧，转动弹簧弹性刚度也取构件弯曲刚度的1000 倍。

图2 角钢构件屈服截面示意图Fig. 2 Schematic diagram of yield section for angle steel member

在进行连续倒塌仿真时，失效构件的变形通常非常大，因此，本文杆件受拉失效准则定义如下：

式中：δ 为构件伸长率；l0和l分别为组合式屈曲单元的初始长度和变形后长度；δT,cr为钢材的临界伸长率，取20%。当构件伸长率超过临界值，认为构件失效，将删除中间的拉压弹簧和弯曲弹簧，此时组合式屈曲单元分成两个独立的弹性梁单元。

该组合式屈曲单元采用拉压弹簧和转动弹簧实现了塑性铰功能，简化了塑性铰法的编程复杂度、降低了刚度矩阵维度、大幅缩短了计算时长。

2 静动力分析方法介绍

2.1 静力分析

静力分析即分析结构在给定静力载荷作用下的响应，并求得结构位移、应力和应变等数据。本文自编求解程序考虑结构的几何非线性，采用更新拉格朗日列式法建立静力平衡方程[20]：

式中：m为单元个数；T为单元从整体坐标系到局部坐标系的转换矩阵； Δu为整体坐标系下的节点位移增量；F为外荷载列阵；为单元弹性刚度矩阵；为单元几何刚度矩阵；为大位移刚度矩阵，是由大位移引起的结构刚度变化，采用增量求解时可忽略[20]。

为避免上述方程隐式迭代计算中由于压杆屈曲失稳带来的刚度突变，进而导致程序不收敛的问题，本文求解静力平衡方程采用显式方法，计算时根据非线性程度将静力平衡计算分为N步，第i步的荷载为Fi=F×i/N，并根据静力平衡方程进行每一步的平衡计算。

2.2 动力分析

在采用本文提出的组合式屈曲单元进行动力分析时，结构整体质量矩阵、刚度矩阵及荷载向量的组装方法均和传统有限元理论一致，不再赘述。采用Rayleigh 公式生成阻尼矩阵时，为避免结构进入非线性后，刚度矩阵变化带来的阻尼矩阵不稳定，本文Rayleigh 阻尼只考虑结构整体质量矩阵。

求解动力平衡方程时采用隐式-显式混合计算方法，隐式算法采用Newmark-β 法，显式算法采用中心差分法。在时程分析前期，无材料非线性行为(未出现塑性铰)，几何非线性影响也较小(位移较小)，隐式算法所需迭代的次数较少，此时，采用隐式算法的计算效率明显高于显式算法。当结构出现强非线性或是杆件发生失稳破坏导致结构整体倒塌时，隐式算法无法收敛，此时，应该采用显式算法。因此，本文编制的求解器首先采用隐式算法进行动力时程分析，当隐式算法不收敛时，将隐式算法最后一个收敛步的结构运动状态作为显式算法初始状态，程序内核调用显式算法继续计算。

3 组合式屈曲单元及求解程序验证

本文提出的组合式屈曲单元需通过自编求解程序进行结构静动力分析，采用MATLAB 程序开发了基于组合式屈曲单元的输电塔结构静动力仿真平台。为说明自编程序的精度及新单元的适用性，本章将与通用有限元分析软件ABAQUS 及实验结果进行全面对比验证。

3.1 数值验证

为验证自编程序在强几何非线性、线弹性状态下，自编程序的刚度矩阵、质量矩阵和有限元求解方法的准确性，本章节将与解析解、通用有限元分析软件进行一系列对比。

3.1.1 解析解对比

用自编程序及ABAQUS 软件分别建立材料、截面、长度均相同的悬臂梁模型，截面为圆环，两个模型杆件均划分10 个单元，且忽略自重影响，具体结构参数详见表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 个单元，并在梁端施加相同大小的集中力计算结构响应，具体结构及荷载参数详见表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 所示，结构底层中心为坐标原点，各节点坐标见表3，同高度各节点x、y坐标绝对值一致，所有构件参数一致，具体见表4。结构前十阶自振频率对比详见表5，最大相对误差只有0.41%，与ABAQUS 软件计算结果非常接近。

表3 框架结构节点坐标Table 3 Node coordinates of frame structure

表4 框架结构构件参数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 号节点x方向位移时程的对比如图6 所示，两条曲线几乎重合，最大误差仅为0.12%，说明自编程序的动力仿真求解器具有较高的计算精度。

图6 框架结构位移时程曲线对比Fig. 6 Comparison of time-history curves of displacement for frame structure

3.2 实验验证

3.2.1 构件滞回特性对比

以Black 等[21]的实验数据为基础，选取4 根钢支撑进行数值模拟，钢支撑具体参数如表6 所示，Brace1 和Brace2 为两端铰接，Brace3 和Brace4为一端刚接、一端铰接。

表6 钢支撑参数[21]Table 6 Parameters of steel brace[21]

图7 为自编程序和实验结果的滞回曲线对比，红色虚线为自编程序计算结果，黑色实线为实验结果。从图中可以看出，自编程序计算结果与实验结果较为吻合，本文提出的组合式屈曲单元能有效地模拟钢支撑受压屈曲后的负刚度现象。

图7 构件滞回曲线对比Fig. 7 Comparison of hysteresis curves of members

图中曲线受拉段的高度不相等，这是由于模拟受拉屈服的拉压弹簧采用理想弹塑性力-位移模型，而实际构件受拉屈服后进入强化阶段，抗拉强度会有所提高。细长构件在外荷载作用下通常发生屈曲失稳，因而，采用理想弹塑性力-位移模型可以很好地模拟该类结构的倒塌破坏行为。

偏心距作为组合式屈曲单元的重要参数，对计算结果有较大影响，现以表6 中Brace1 杆件为例，采用不同的偏心距参数计算最大受压承载力，详见表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 所示，图中圆点表示塑性铰位置，塑性铰均出现于塔腿处，与实验塔破坏位置相同。自编程序模拟的承载能力与实验结果非常接近，进一步验证了本文提出的组合式屈曲单元的准确性和适用性。对比实验和仿真的连续倒塌动画和倒塌破坏位置可知，各时刻的倒塌状态非常相近，说明了本文开发的动力仿真平台在连续倒塌模拟方面具有很高的精度。

图9 推覆分析中输电塔塑性铰位置Fig. 9 Position of plastic hinges of transmission tower in pushover analysis

为说明本文提出的组合式屈曲单元及自编程序在计算效率方面的优势，现以上述输电塔结构为例，自编程序生成的刚度矩阵维度为24 519×24 519，通用有限元软件为了保证计算精度，将每个杆件划分为5 个单元，生成刚度矩阵维度为383 721×383 721。同时以100%的参考荷载为例，自编程序计算时长3.3254 s，通用有限元软件计算时长25.501 25 s，电脑CPU 型号为i5-8500，内存16 GB。由此可见，本文提出的组合式屈曲单元及自编程序可大幅度提高计算效率，且具备较高计算精度。

4 结论

本文基于有限元理论，提出了一种组合式屈曲单元模拟细长杆件的屈曲失稳现象，囊括了三种常见的约束类型：

(1)该单元采用拉压弹簧模拟单元轴向塑性伸长，采用转动弹簧模拟单元塑性弯曲，实现了塑性铰功能，简化了塑性铰法的编程复杂度，降低了刚度矩阵维度，大幅缩短了计算时长。

(2)基于MATLAB 程序编写了该单元的静动力仿真程序，可以开展工程结构的静力推覆分析和动态倒塌仿真分析。

(3)通过和解析解、ABAQUS 商用软件、构件实验、输电塔足尺实验对比，全面验证了该组合式屈曲单元的准确性及自编仿真程序的可靠性。