单层柱面网壳动力稳定性评估的阶跃推覆分析方法

曲 扬 罗永峰 黄青隆

(同济大学土木工程学院，上海 200092)

网壳结构的动力失稳是一个非线性动力学问题，是结构地震反应分析和抗震设计的关键问题之一.网壳结构的失稳是结构由高位能向低位能快速转变的过程，伴随着结构位形的动态变化，使得结构的失稳具有动态特征.因此，网壳结构的动力稳定性研究具有广泛意义[1].

对网壳结构稳定性的研究始于静力稳定问题.目前，特征节点的荷载-位移曲线是最常用的静力稳定分析曲线，但该曲线仅能反映单个节点的受力变形路径，且选取特征节点主观性强，对于体系复杂、受力变形关系并不直观的网壳结构而言，该曲线的缺点愈发明显[2].为全面掌握网壳结构整体稳定性，学者们提出了切线刚度最小特征值[3]、整体刚度参数[4]、特征刚度[2]等性能参数，用以反映结构在受载过程中刚度的变化，能够较好地克服特征节点荷载-位移曲线的缺点，并得以广泛应用和推广[5-6].相比静力稳定问题，网壳结构的动力稳定特性更加复杂，若仍然沿用静力稳定分析方法判定动力失稳将可能导致误判.文献[7]指出，拟静力刚度准则不足以判定动力失稳，往往会低估结构动力失稳临界荷载，建议采用能量准则判定结构动力失稳.在此基础上，文献[8-10]从能量变化的角度出发，研究了单层网壳结构的动力稳定性，并提出了相应的定量判定准则.文献[11-12]基于结构响应突变准则，对单层柱面网壳的强震失效机理进行了研究，归纳了动力失稳和强度破坏2种失效模式，并进一步提出了模糊综合判定方法[13].上述判定准则虽然可以从能量角度给出动力失稳的合理解释，但所采用的时程分析方法需要大量计算才能获得网壳结构动力失稳临界荷载，实用性不强.近年来应用广泛的推覆分析方法虽可用于评估网壳结构地震响应[4]，但由于推覆曲线后屈服段较短[14]，且无法全面考虑动力效应，因此其应用范围仅限于强度破坏模式，难以直接用于评估网壳结构动力稳定性.

针对上述问题，本文采用推覆分析构造基于振型刚度的ESDF体系，确定非线性能力曲线，然后对网壳结构进行逐级阶跃加载推覆分析，根据其临界状态荷载因子，对ESDF体系进行阶跃加载时程分析，结合地震反应谱，最终获得网壳结构动力失稳临界加速度因子，从而建立阶跃推覆分析方法.将此方法运用于某单层柱面网壳算例的动力稳定性评估，并将计算结果与时程分析方法和传统推覆方法结果进行对比，以验证本文方法的适用性与计算效率.

1 阶跃推覆分析方法

1.1 基于振型刚度的ESDF体系

是否能够构造合理的等效单自由度(ESDF)体系，决定了能否准确求解网壳结构的地震响应.传统推覆分析方法采用基底剪力和顶点位移构造结构的能力曲线，无法直观准确地反映网壳结构的受力变形关系[15].鉴于此，本文采用基于振型刚度的改进模态推覆分析方法[16]，构造主振型ESDF体系.

对于具有N个节点的网壳结构，在地震作用üg(t)激励下，其非线性动力微分方程为

(1)

式中，u(t)为地震作用下的结构动力位移；M为质量矩阵；C为阻尼矩阵，采用经典阻尼；ι为影响因子向量，地震激励所在方向上的元素取1，其他方向上的元素取0；Fs(t)为非线性恢复力，与当前变形状态和加载历史相关.

将线性体系中的模态分解方法引入非线性分析中，可将非线性体系的动力位移分解为

(2)

式中，φl，ql(t)分别为第l阶振型的振型向量和坐标.

(3)

式中，ωn，ζn分别为第n阶振型频率和阻尼比；mn，Ln分别为振型质量和参量，可表示为

(4)

对于第n阶振型，qn(t)可以写为qn(t)=ΓnDn(t)，其中Dn(t)为振型位移.将qn(t)=ΓnDn(t)代入式(3)，并考虑到振型参与系数Γn=Ln/mn，第n阶振型动力方程可解耦成以Dn(t)为基本未知量的微分方程，即

(5)

令伪加速度为

(6)

(7)

根据推覆分析原理，对于非线性体系，第i荷载步的推覆荷载增量ΔPi和位移响应增量Δui关系由增量平衡方程给出，即

ΔPi=KT,iΔui

(8)

式中，KT,i为第i荷载步的结构切线刚度矩阵.

非线性体系第n阶振型刚度定义式为[16]

(9)

由式(9)可知，振型刚度能够反映结构整体刚度矩阵的变化，且与刚度矩阵量纲一致，易于理解，便于计算.此外，采用振型刚度描述结构非线性特性能够考虑结构整体响应，避免片面地依赖于特定节点和特征响应.

通过理论推导可知，非线性体系第n阶振型的刚度与本文采用的振型刚度在数值上近似相等，且量纲一致.因此，伪加速度增量ΔAn,i和振型位移增量ΔDn,i可表示为

(10)

采用基于振型刚度构造的ESDF体系，求解式(7)可得到振型响应.将各阶振型响应组合，最终得到结构总响应.

1.2 阶跃推覆分析与失稳临界状态

为在推覆分析方法中引入动力效应，可对网壳结构逐级施加阶跃荷载，进行推覆分析.对于有阻尼非线性多自由度体系，按照第n阶振型模式施加阶跃荷载χMφn，其动力方程为

(11)

式中，χ为荷载因子.

在任意激励作用下，结构动力稳定性判别可采用Budiansky-Roth准则[17].该准则通过在逐级加载计算中观察结构体系的动力响应变化，判定体系的动力稳定性.将增量计算过程中荷载的微小增量可能导致结构响应突然增大时的结构与荷载状态定义为失稳临界状态，此时对应的荷载定义为结构动力失稳临界荷载.

根据Budiansky-Roth准则，可控制荷载因子χ对结构进行逐级阶跃加载推覆分析，从而得到结构在阶跃荷载χMφn作用下第n阶临界状态荷载因子χn.对网壳施加阶跃荷载χnMφn时，结构达到失稳临界状态，若继续增大荷载，结构将发生动力失稳.

为考察网壳结构在阶跃荷载作用下的第n阶振型响应，将线性体系中的模态分解原理引入非线性分析中，则式(11)解耦为

(12)

式中，Γn为第n阶振型参与系数.

式(12)为第n阶振型ESDF体系在临界阶跃荷载χn/Γn作用下的基本动力方程.由于网壳结构在动力失稳之前，整体刚度衰减较小，仍能保持基本完整的受力形态[18]，因此，ESDF体系仍采用1.1节中的非线性能力曲线.采用时程分析方法求解式(12)，即可得到ESDF体系在临界阶跃荷载χn/Γn作用下的位移时程Dn(t)，进而得到伪加速度峰值An,max.网壳结构在临界阶跃荷载χnMφn作用下，第n阶振型ESDF体系的伪加速度峰值为An,max.

1.3 动力失稳临界加速度因子

为衡量网壳结构在不同地震动强度下的动力响应，需对地震动输入进行调幅.现定义加速度因子λ为调幅系数，则结构在λüg(t)作用下，第n阶振型伪加速度峰值为λSa(Tn)，其中Sa(Tn)为谱加速度值.假定当λ=λc时，网壳结构达到动力稳定临界状态，若继续增大地震动强度，结构将发生动力失稳，则λc为动力失稳临界加速度因子.

由前述可知，网壳结构在阶跃荷载作用下，临界状态的第n阶伪加速度峰值为An,max；在地震作用下，临界状态的第n阶伪加速度峰值为λSa(Tn).因此，基于网壳结构在临界阶跃荷载作用下的失稳临界状态，评估在地震作用下的动力稳定性，需满足

An,max=λnSa(Tn)

(13)

式中，λn为假定网壳结构按第n阶振型模式失稳时的临界加速度因子.

对于各阶主振型，计算相应的λn可得到临界加速度因子的集合{λn}，取最终的动力失稳临界加速度因子λc=min{λn}.

1.4 计算步骤

本文提出的阶跃推覆方法计算步骤如下：

① 对网壳结构进行模态分析，获得各阶振型频率ω和振型参与系数Γ，选取主振型；

② 针对选取的第n阶主振型进行模态推覆分析，构造相应的ESDF体系，获得非线性能力曲线；

③ 对网壳结构进行第n阶阶跃推覆分析，计算得到临界状态荷载因子χn；

④ 对第n阶ESDF体系施加阶跃荷载χn/Γn，并进行时程分析(即求解式(12))，得到对应的临界失稳伪加速度峰值An,max；

⑤ 针对选取的各阶主振型，重复步骤②～步骤④，得到对应各阶主振型的临界失稳伪加速度；

⑥ 结合地震反应谱，根据式(13)，求得对应各阶振型的失稳临界加速度因子λn，取最终的动力失稳临界加速度因子λc=min{λn}.

2 柱面网壳算例分析

为验证本文提出的阶跃推覆分析法的适用性及其计算效率，采用ANSYS和MATLAB对某一典型柱面网壳结构算例进行静力推覆分析、阶跃推覆分析和时程分析.以不同方法失效时刻的结构位移响应、单元应力响应作为主要参数，考察阶跃推覆分析法的合理性，并对比各方法失效时刻的临界荷载，分析阶跃推覆分析法的误差和效率.

2.1 结构模型参数

采用三向网格型单层柱面网壳结构作为算例(见图1).该网壳跨度为15 m，矢跨比为1/2，长宽比为1.0，采用焊接球节点，纵边支承为固定铰支座，端部山墙支承为竖向铰支座.网壳屋面均布荷载为1.5 kN/m2，将其转化为等效集中质量施加于节点，以考虑其惯性力效应.在满足静力设计要求下，结构构件均采用圆钢管，横杆和纵杆均为φ89 mm×4 mm，斜杆φ114 mm×4 mm.构件材料为Q235钢，采用双线性随动强化准则，弹性模量为206 GPa，屈服强度为235 MPa，屈服后弹性模量为0.8 GPa.采用Rayleigh阻尼，阻尼比为0.02.采用通用有限元软件ANSYS中的Beam189单元和Mass21单元模拟结构构件和节点集中质量，计算分析中计入结构的材料非线性和几何非线性效应.

2.2 主振型ESDF体系

为快速准确选取网壳结构主振型，本文采用Luo等[19-20]提出的振型遴选阈值法识别主振型.基于静力功参与系数(SPF)和动力功参与系数(DPF)，采用参考阈值δ=0.20作为振型遴选标准，即选择参与系数SPF和DPF均大于0.20的振型作为主振型.表1和图2分别给出了基于阈值法遴选得到的主振型基本信息和振型图.对网壳结构按照主振型进行非线性模态推覆分析，得到ESDF体系振型刚度k*、伪加速度A和振型位移D的变化，进而绘制出非线性能力曲线.图3给出了主振型ESDF体系的k*-D曲线和A-D曲线.

表1 网壳结构主振型基本信息

图2 网壳结构竖向主振型图

(a) k*-D曲线

(b) A-D曲线

2.3 Taft地震波计算结果对比

为初步验证本文方法预测动力失稳临界点的有效性，同时全面考察该方法给出的动力失稳模式的合理性，采用1952年美国KernCounty地震中记录的Taft地震波竖向分量作为地震动输入，详细对比时程分析方法、阶跃推覆分析方法和传统推覆分析方法计算得到的网壳结构动力失稳临界点及其地震响应.

为考察单层柱面网壳结构在Taft地震波作用下的动力失稳临界点，将3种方法计算得出的伪加速度An与最大位移Dmax关系曲线绘于图4，并在图中标出相应的动力失稳临界点.其中，传统推覆方法以最后一步计算的收敛点作为失稳临界点.计算结果表明，时程分析方法和阶跃推覆方法计算得到的最大位移值在发生动力失稳之前随伪加速度值变化缓慢，网壳结构仍有较大刚度；当伪加速度达到约3.5g时，结构最大位移仅为500～700 mm；当伪加速度在此水平上稍加增大，结构位移响应即发生明显的突然增大现象，属于典型的动力失稳破坏.时程分析方法与阶跃推覆方法计算得到的动力失稳临界点较为接近，说明阶跃推覆分析方法由于考虑了动力作用，可以较为准确地预测网壳结构动力失稳临界荷载.相比之下，纯静力的传统推覆方法由于推覆分析后段不易收敛，造成推覆曲线后屈服段较短，且其An-Dmax曲线不具备动力失稳破坏特征，失效临界荷载显著小于真实的结构动力失稳临界荷载，说明传统推覆方法无法模拟出动力失稳破坏现象，难以评估网壳结构动力稳定性.

图4 An-Dmax曲线对比

为考察3种方法给出的动力失稳临界点处的动力响应，网壳结构的竖向位移响应uz和单元最大应力响应σmax对比见图5和图6.其中，单元最大应力响应取单元各积分点的应力最大值.由图可知，时程分析方法和阶跃推覆方法在动力失稳临界点处的变形模式和单元应力分布规律较为一致，较大挠度和高应力区均集中于中部跨中区域，且发展程度接近，表明阶跃推覆方法可以在一定程度上预测网壳结构的动力失稳模式.对应的传统推覆方法在动力失稳临界点处的结构变形和单元应力响应均显著小于时程分析方法的计算结果，说明纯静力的传统推覆方法低估了网壳结构在动力失稳前后的位移响应和单元应力响应发展程度.

(a) 时程分析方法

(b) 阶跃推覆方法

(c) 传统推覆方法

(a) 时程分析方法

(b) 阶跃推覆方法

(c) 传统推覆方法

2.4 不同地震波计算结果对比

为进一步验证阶跃推覆分析方法对不同地震动波的适用性，本文从日本防灾科学技术研究所提供的K-NET强震数据库中选取了近20年来发生的13次6.0级以上地震中共计45条地震动记录的竖向分量(E1～E45)，作为结构的地震动输入，其伪加速度反应谱见图7，详细信息见表2.

表2 地震波详细信息

图7 地震波伪加速度反应谱

采用时程分析方法、阶跃推覆方法以及传统推覆方法计算网壳结构的失稳临界加速度因子λc，并绘制于图8中.分析结果表明，传统推覆方法明显低估了网壳结构的动力失稳承载力，对比时程分析方法，平均误差达56.7%，计算结果过于保守且精度较低，无法准确预测网壳结构动力失稳临界加速度因子.阶跃推覆方法计算结果与时程分析方法计算结果的误差基本控制在40%以内，平均误差仅为23.6%，相比传统推覆方法，计算精度提高了58.3%，可以较为快速准确地评估网壳结构动力稳定性.

图8 失稳临界加速度因子对比

2.5 计算效率对比

与时程分析方法相比，阶跃推覆方法能够保证计算精度，同时大幅提高计算效率，分析耗时仅约为时程分析方法的25%.与传统推覆方法相比，阶跃推覆方法尽管分析耗时增加了约80%，但计算精度提高了58.3%.若结合合理的弹塑性反应谱，阶跃推覆方法的计算精度将进一步提高.因此，本文提出的阶跃推覆方法能够以较低的计算代价，得到满足工程需求的计算精度，可用于快速准确评估单层柱面网壳结构在地震作用下的动力稳定性.

3 结论

1) 与纯静力的传统推覆分析方法相比，本文方法可以考虑动力效应，能够有效模拟单层柱面网壳结构的动力失稳现象.

2) 本文方法计算得到的动力失稳临界点处的变形模式和受力状态与时程分析方法结果较为一致，基本能够预测单层柱面网壳结构的动力失稳模式.

3) 本文方法能够较为准确地预测单层柱面网壳结构的动力失稳临界加速度因子，平均误差仅为23.6%，计算代价小，分析效率高，是一种快速准确评估单层柱面网壳动力稳定性的实用方法.