0 引言

　　 空间桁架体系不仅被大量运用到大型体育馆、机场、桥梁等大跨结构中,而且在输电塔、电视塔等高耸结构中也应用广泛。该类型的结构体系传力路径明确、受力合理,但对材料非线性、几何非线性敏感。往往单根空间杆件的失效就会引起整体结构的动力响应,甚至造成结构整体倒塌。因此,对强烈地震作用下桁架结构的地震倒塌进行计算和分析具有重要意义。

　　 近年来,国内外学者均对空间桁架结构的倒塌进行了一些研究。杨明飞、徐赵东 ^[1]基于有限元软件,同时考虑结构的双重非线性和材料失效应变,对三种典型大跨网格的倒塌过程进行了模拟,并提出了大跨网格结构的倒塌模式。朱奕锋等 ^[2]采用去除杆件法对某会展中心进行了连续倒塌分析,确定了该结构承载的关键部位。江晓峰、陈以一 ^[3]通过LS-DYNA模拟了一重型钢结构厂房屋架的连续倒塌过程,归纳了三种基本的内力重分布机制。赵广坡 ^[4]等对某大跨屋架结构进行了静力非线性分析,并论述了该结构的抗倒塌设计方法。

　　 在地震作用下,桁架结构中的杆件受压易发生屈曲; 反向作用时杆件又会突然卸载,甚至受拉屈服,这些过程均会引起结构内力重分布而产生动力效应 ^[5,6]。因此,空间桁架结构的地震响应分析需要考虑杆件的这些非线性行为 ^[7]。谢道清等 ^[8]拟合了一个同时考虑受拉屈服和受压屈曲的圆管杆单元模型,并把该模型运用于一球面网壳结构的罕遇地震时程计算。Thai等 ^[9]采用Hill杆件屈曲模型 ^[10],模拟了强震作用下某空间桁架结构的动力响应,并与理想弹塑性模型的结果进行了对比,证明了杆件屈曲模型开发的必要性。Zheng等 ^[11]采用不同杆单元模型评估了某桁架结构抗地震倒塌能力,结果表明可模拟压屈的杆单元模型可以更加合理地评估结构的抗倒塌能力。

　　 从上述可知,有限元软件ANSYS和ABAQUS是常用来模拟桁架结构倒塌的工具,但这些软件中没有能考虑受压屈曲行为的杆单元,因此不能准确模拟构件受压屈曲后结构的真实响应。一些学者虽然编制了压屈单元程序,但没有考虑杆件断裂破坏对结构的影响。因此,本文首先提出了一种考虑杆件受拉屈服、受压屈曲行为的杆件压屈模型; 然后定义了杆件断裂破坏准则,并嵌入到显式动力计算程序中,实现了空间桁架结构的地震倒塌模拟。

　　 图1 强度破坏准则   

　　  

　　 图2 极限应变破坏准则  

　　  

　　 图3 压杆屈曲破坏准则  

　　  

1 杆件破坏准则

1.1 强度破坏准则

　　 对脆性材料,杆件的断裂破坏可采用极限应力作为指标,如图1所示。当杆件的拉压应力在达到材料的极限应力之前,加载、卸载均为线弹性; 当杆件拉压应力达到极限应力后,杆件突然断裂,应力变为零。杆件的应力-应变关系为:

　　 σ={Eε0(−εcr<ε<εy)(ε≤−εcr或ε≥εy)         (1)σ={Eε(-εcr<ε<εy)0(ε≤-εcr或ε≥εy)         (1)

　　 其中:

　　 ε_y=σ_y/E; ε_cr=min(σ_y/E,σ_cr/E)

　　 式中:σ_cr=π²E/λ²为杆件的欧拉临界应力(图1), λ为杆件的长细比; ε为杆单元的轴向应变; E为杆件的弹性模量; σ_y为杆件的屈服应力。

1.2 极限应变破坏准则

　　 对延性较好、长细比较小的杆件,其断裂破坏可采用极限应变作为指标。这类材料的应力-应变模型可简化为双线性模型,如图2所示。这类模型拉压对称,弹性卸载,当拉压应变达到极限应变时杆件发生断裂,其应力-应变关系为:

　　 σ=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪Eε        (−εcr<ε≤εy)σy+Et(ε−εy)  (εy<ε≤εcrit)−σy+Et(ε+εy) (−εcrit<ε≤−εy)0        (ε≤−εcrit或ε>εcrit)(2)σ={Eε        (-εcr<ε≤εy)σy+Et(ε-εy)  (εy<ε≤εcrit)-σy+Et(ε+εy) (-εcrit<ε≤-εy)0        (ε≤-εcrit或ε>εcrit)(2)

　　 式中:ε_crit为杆件轴向断裂的极限应变;E_t为杆件的塑性模量。

1.3 压杆屈曲软化破坏模型

　　 对延性较好、长细比较大杆件,其受压时会发生非线性屈曲,导致受力杆件突然丧失承载力,引起结构动力响应。因此,若要准确模拟结构动力响应,杆件单元必须能够考虑受压非线性屈曲行为。本文对Hill模型 ^[10]进行了改进,改进后的模型不仅可以模拟杆件在往复荷载作用下的加载、卸载规律,还可以考虑杆件受拉屈服和受压非线性屈曲行为,见图3。

1.3.1 加载阶段

　　 若杆件初始受拉,则应力-应变关系为:

　　 σ=⎧⎩⎨⎪⎪Eε  (0<ε≤εy)σy  (εy<ε≤εcrit)0  (ε>εcrit)         (3a)σ={Eε  (0<ε≤εy)σy  (εy<ε≤εcrit)0  (ε>εcrit)         (3a)

　　 若杆件初始受压,则应力-应变关系为:

　　 σ={Eε  (−εcr<ε≤0)σcr  (−εb<ε≤−εcr)         (3b)σ={Eε  (-εcr<ε≤0)σcr  (-εb<ε≤-εcr)         (3b)

　　 式中-ε_b为B点所对应的应变。

　　 大量试验表明,非线性屈曲软化BC段(图 3)的行为与杆件的截面形式及长细比相关。Hill模型 ^[10]采用指数函数描述该非线性屈曲段,但函数中的参数X₁和X₂被指定为两个常数,未考虑杆件截面形式及长细比的影响。因此,本文基于塑性铰理论对该段进行了推导。大量试验表明当杆件受压发生非线性屈曲时,跨中处已形成塑性铰 ^[12],如图4(a)所示。跨中截面需满足轴力-弯矩屈服面函数,本文屈服面采用AISC-LRFD双线性关系式:

　　 ⎧⎩⎨⎪⎪P2Py+MprMp=1  (PPy<0.2)PPy+8Mpr9Mp=1  (PPy≥0.2)         (4){Ρ2Ρy+ΜprΜp=1  (ΡΡy<0.2)ΡΡy+8Μpr9Μp=1  (ΡΡy≥0.2)         (4)

　　 式中:P为轴向力; M_pr为跨中弯矩; M_p为P=0时对应的极限塑性弯矩; P_y为屈服荷载。

　　 根据弯矩平衡方程可得:

　　 Mpr=PΔ         (5)Μpr=ΡΔ         (5)

　　 式中Δ为跨中处的横向位移,见图 4(b)。

　　 杆件轴向总位移δ主要由两部分组成:弹性变形引起的轴向位移δ_cr和杆件弯曲变形(几何非线性)引起的轴向位移δ_b,即:

　　 δ=δcr+δb         (6)δ=δcr+δb         (6)

　　 式中δ_cr=σ_crL/E,其中L为杆件长度。

　　 杆件弯曲变形可做如图 4(c)简化。根据简化的变形图,利用勾股定理可得:

　　 Δ=(L2)2−(L2−δb2)2−−−−−−−−−−−−−−−√         (7)Δ=(L2)2-(L2-δb2)2         (7)

　　 对式(7)整理后可得杆件弯曲引起的轴向变形δ_b为:

　　 δb=L−L2−4Δ2−−−−−−−−√         (8)δb=L-L2-4Δ2         (8)

　　 联立式(4)～(8),可得:

　　 σ=⎧⎩⎨⎪⎪⎪⎪2MPσyPyL1−(ε−εcr+1)2√+Mp  (σσy<0.2)9MPσy4PyL1−(ε−εcr+1)2√+9Mp (σσy≥0.2)         (9)σ={2ΜΡσyΡyL1-(ε-εcr+1)2+Μp  (σσy<0.2)9ΜΡσy4ΡyL1-(ε-εcr+1)2+9Μp (σσy≥0.2)         (9)

　　 式中:σ=P/A; ε=δ/L; ε=δ_cr/L。

　　 式(9)即为BC段对应的应力-应变关系。由于B点处(图 3)的应力与应变也满足式(9),且已知B点处的应力等于σ_cr,代入式(9)得:

　　 εb=⎧⎩⎨⎪⎪⎪⎪⎪⎪1−1−[Mp(2σy/σE−1)PyL]2−−−−−−−−−−−−−−−√+εcr  (σσy<0.2)1−1−[9Mp(σy/σE−1)4PyL]2−−−−−−−−−−−−−−−√+εcr  (σσy≥0.2)         (10)εb={1-1-[Μp(2σy/σE-1)ΡyL]2+εcr  (σσy<0.2)1-1-[9Μp(σy/σE-1)4ΡyL]2+εcr  (σσy≥0.2)         (10)

　　 对式(10)整理后可得:

　　 εb−εcr=⎧⎩⎨⎪⎪⎪⎪⎪⎪1−1−[Mp(2σy/σE−1)PyL]2−−−−−−−−−−−−−−−√  (σσy<0.2)1−1−[9Mp(σy/σE−1)4PyL]2−−−−−−−−−−−−−−−√  (σσy≥0.2)         (11)εb-εcr={1-1-[Μp(2σy/σE-1)ΡyL]2  (σσy<0.2)1-1-[9Μp(σy/σE-1)4ΡyL]2  (σσy≥0.2)         (11)

　　 式中ε_b-ε_cr的值为AB平台段长度的大小(图 3)。

1.3.2 卸载-再加载阶段

　　 若杆件处于受拉状态,卸载-再加载路径均为弹性过程; 若杆件处于受压状态,卸载-再加载路径取决于当前的应力水平。当在OA,AB段卸载时,其为弹性过程; 当在BC段卸载时,卸载路径非常复杂,很难得到一个精确的本构方程,本文假设应力-应变沿直线返回到屈服应力的一半处 ^[10],如图 3中CD段所示,其应力-应变表达式为:

　　 σ=σD−σCεD−εC(ε−εC)+σC         (12)σ=σD-σCεD-εC(ε-εC)+σC         (12)

　　 式中:σ_D,ε_D分别为D点的应力和应变; σ_C,ε_C分别为卸载起点C的应力和应变。

　　 图4 杆件压屈后变形图  

　　  

2 显式动力分析方法

　　 三维空间桁架结构(图 5(a))可以被离散为一系列质点和杆元(图 5(b)),结构的质量平均分配到每个质点上。结构在运动过程中,变形后的杆件可以简化为一系列力作用在质点上,如图5(c)所示。

　　 图5 空间桁架简化模型  

　　  

2.1 运动方程及求解

　　 在地震动下,某质点i的运动方程可表示为:

　　 mid⋅⋅=Fexti−miαg−Finti−miξd˙         (13)mid⋅⋅=Fiext-miαg-Fiint-miξd˙         (13)

　　 式中:d⋅⋅,d˙d⋅⋅,d˙分别为加速度向量和速度向量; Fextiiext,Fintiiint分别为外力向量和内力向量,其中Fintiiint可从与质点i相连的杆件获得; m_i为质点i的质量; ξ为阻尼系数; α_g为三维地震地面运动加速度向量。

　　 为避免大量迭代计算和方便处理断裂问题,式(13)采用中心差分法求解:

　　 d˙=12Δt(dn+1−dn−1)         (14a)d⋅⋅=1Δt2(dn+1−2dn+dn−1)         (14b)d˙=12Δt(dn+1-dn-1)         (14a)d⋅⋅=1Δt2(dn+1-2dn+dn-1)         (14b)

　　 式中:d_n+1,d_n,d_n-1分别为n+1,n,n-1时刻杆件的轴向位移; Δt为时间步长,本文取1×10^-4s。

　　 把式(14)代入式(13)后可得第n+1时间步质点i的位移为:

　　 dn+1=2C1dn−C2dn−1+C1Δt2mi(Fexti,n−Finti,n−miαg,n)         (15)dn+1=2C1dn-C2dn-1+C1Δt2mi(Fi,next-Fi,nint-miαg,n)         (15)

　　 式中:C1=1/(1+12ξΔt)C1=1/(1+12ξΔt);C2=1/(1−12ξΔt)C2=1/(1-12ξΔt)。

2.2 杆件断裂后处理

　　 在动力分析过程中,当杆件达到图1～3所定义的断裂准则时,假设杆件在两端发生破坏,如图6所示。∣∣FintA∣∣|FAint|和∣∣FintB∣∣|FBint|为A,B两点内力的大小,断裂模式如下:当∣∣FintA∣∣>∣∣FintB∣∣|FAint|>|FBint|时,杆件在A端断裂,如图6(b)所示; 当∣∣FintA∣∣<∣∣FintB∣∣|FAint|<|FBint|时,杆件在B端断裂,如图6(c)所示; 当∣∣FintA∣∣=∣∣FintB∣∣|FAint|=|FBint|时,杆件在A,B端同时断裂,如图6(d)所示。

3 数值算例

3.1 算例1

　　 在该三层空间桁架中,竖向杆件、水平杆件及斜撑均为圆形钢管,其质量分布、几何和材料性质及桁架模型如图 7所示,每个结点处的集中质量为10 000kg(实心点处)。

　　 图6 断裂模型  

　　  

　　 图7 空间桁架模型  

　　  

　　 为方便与有限元软件的结果进行比较,本文基于理想弹塑性杆单元模型(Ideal Elastic-plastic Model, IEM),分别采用ANSYS和显式动力分析方法计算了该三层空间桁架结构在El Centro波作用下的动力响应。图8比较了结构顶点处的时程位移曲线,可发现显式算法与ANSYS的计算结果吻合较好,响应位移的最大误差为5.02%; 图9比较了底层杆件1(图 7)在地震动作用下的滞回曲线,除在曲线拐角处存在偏差,其余地方两者基本一致,以上对比证明本文自编程序具有可靠性。

　　 图8 理想弹塑性模型下位移时程曲线  

　　  

　　 图9 滞回曲线(理想弹塑性单元)  

　　  

　　 通用有限元软件缺少可以模拟杆件受压非线性屈曲的杆单元模型(Nonlinear Buckling Member Model, NIBM), 因此本文将压杆屈曲单元(图 3)嵌入到显式动力程序中,模拟了结构的地震动力响应,并与采用理想弹塑性单元模型的计算结果进行了对比,如图10所示。图11对比了杆件1的滞回曲线。结果表明,结构中的空间杆件在地震作用下发生了明显的非线性屈曲,结构顶点响应位移也较大。因此,在强震作用下,采用理想弹塑形模型已不能准确预测该结构的动力响应,同时也证明了压杆屈曲模型开发的必要性。

　　 图10 不同单元模型下位移时程曲线 

　　  

　　 图11 不同单元模型下滞回曲线   

　　  

3.2 算例2

　　 某高耸输电塔总高122m,呼称高110m,根开20.24m,塔架截面为正方形,结构总质量约为1.3×10⁵kg。输电塔主材、横隔和斜撑采用薄壁钢管,该塔有830个杆单元和253个质点。输电塔主材、辅材、斜撑均采用Q235钢材,材料的屈服强度为235MPa,弹性模量200GPa,质量密度为7 800kg/m³,泊松比为0.3。输电塔-线体系由输电塔、导线和绝缘子组成,与输电塔的质量相比,导线的质量较小,为简化计算,导线简化为等效质量,将其添加到与输电塔连接的部位,且每个集中质量为1×10⁴kg。输电塔模型及等效质量分布如图 12 所示。

　　 图12 输电塔及导线简化模型  

　　  

　　 若杆件采用强度破坏准则(图 1),结构发生倒塌的最小加速度峰值PGA为9.8m/s²。图13(a)给出了结构倒塌过程中塔顶水平和竖向位移时程曲线以及特殊时刻结构的变形图。当t=4.9s时,塔头部分杆件受压达到极限应力首先发生局部破坏; 当t=6.1s时,塔高102m处主材受压发生较大变形,随后此处杆件开始发生断裂,输电塔的水平和竖向位移开始逐渐增大; 当t=7.4s时,塔高102～106m处杆件发生局部破坏; 当t=10.9s时塔头和塔身都有很大程度的破坏,输电塔不能再承担水平和竖向荷载。

　　 图13 塔顶位移时程曲线及特殊时刻变形图 

　　  

　　 若采用压杆屈曲软化破坏准则(图 3,取ε_crit=0.003),结构发生连续倒塌的最小加速度峰值PGA为11.7m/s²,图13(b)给出了结构倒塌过程中塔顶水平和竖向位移时程曲线以及特殊时刻结构的变形图。当t=5.5s时,塔头部分杆件受压屈曲首先发生局部破坏,由于结构杆件的塑性得到了适当发挥,因此破坏时间较强度破坏准则晚; 当t=6.4s时,塔腿发生破坏,输电塔的水平和竖向变形开始逐渐增大; 当t=7.5s时,塔高0～66.5m处较多杆件发生屈曲,输电塔丧失水平和竖向荷载传递路径,塔身上部和塔头在重力作用下坠落到地面。

　　 若杆件采用极限应变破坏准则(图 2),结构发生连续倒塌的最小加速度峰值PGA为14.6m/s²,图13(c)给出了结构倒塌过程中塔顶水平和竖向位移时程曲线以及特殊时刻结构的变形图。当t=6.2s时,塔高22.25m处的杆件受压达到极限应变,首先发生断裂; 当t=8.0s时,塔高22.25m处局部发生较大变形,输电塔的水平和竖向变形开始逐渐增大; 当t=10.2s时,塔高22.25～66.5m处较多杆件达到极限应变值,输电塔丧失水平和竖向荷载传递路径,塔身上部和塔头部分坠落到地面。

　　 通过上述分析可知,采用不同的杆件破坏准则对输电塔结构的地震倒塌计算结果有很大影响。输电塔结构中存在大量长细比较大的支撑杆件,该类型杆件在受压达到屈曲荷载时发生非线性屈曲,并未断裂,因此强度破坏准则低估了结构的抗地震倒塌能力,且结构破坏突然; 而极限应变破坏准则下,结构从局部破坏到整体结构倒塌有一定时间间隔,结构的整体塑性能力得到发挥,但不能考虑长细比较大杆件的非线性屈曲行为,因此将高估结构的抗地震倒塌能力。当采用压杆屈曲软化破坏准则时,结构倒塌的最小加速度峰值在这两种破坏准则之间,该准则既能模拟杆件受拉屈服,又可以考虑受压非线性屈曲,因此可以更为合理地模拟结构的地震倒塌过程。

4 结论

　　 本文提出了一种较精确的杆单元模型,并嵌入到显式动力程序中,编制了空间桁架结构的地震倒塌计算程序。通过对一个三层空间桁架和输电塔结构进行弹塑性分析和地震倒塌分析,得到以下结论:

　　 (1)所提出的压杆屈曲单元模型,既可以模拟杆件受拉屈服,又可以考虑杆件受压非线性屈曲。当分别采用压杆屈曲模型和理想弹塑性单元模型对结构进行非线性分析时,结构响应差别较大,证明了压杆屈曲软化模型开发的必要性。

　　 (2)显式分析方法在计算过程中,不需要形成总体刚度矩阵和增量形式的切线矩阵。当杆件发生断裂时,只需增加相应的质点和方程个数即可,因此能较方便地处理断裂问题。

　　 (3)采用不同的杆件破坏准则将影响输电塔的地震倒塌特征,采用强度破坏准则和极限应变破坏准则将分别低估和高估输电塔的抗震能力,压杆屈曲软化破坏准则既能模拟压杆的屈曲软化破坏,又能模拟拉杆的极限应变破坏。因此,采用该杆件破坏准则能较合理地评价输电塔结构的抗地震倒塌能力。