对结构进行全过程倒塌仿真可有效发现结构发生破坏时的薄弱部位。现有倒塌仿真分析有4种方法, 即:1) 修正有限元法, 主要用于静力分析;2) 离散单元法, 不要求满足连续性条件;3) 直接采用有限元分析软件进行倒塌仿真, 如LS-DYNA;4) 在通用有限元软件 (如MARC, ANSYS, ABAQUS) 平台上进行二次开发。由于房屋建筑结构的体量巨大, 单元之间以共享节点的方式满足连续条件, 计算方式为动力时程分析, 因此前两种方法显然不适用。倒塌分析涉及到各种类型构件失效问题研究, 第4种方法允许用户根据工程需要自行编制计算程序源代码, 以精确控制单元失效, 故相对于第3种方法能解决更为复杂的实际工程难题。

　　 陆新征等^[1]在MARC平台上进行二次开发, 对多个建筑结构进行倒塌仿真, 但是MARC软件采用隐式算法, 计算耗时随模型自由度呈指数增加。房屋建筑结构模型的自由度数量巨大, ABAQUS软件显式算法求解器计算耗时随模型自由度呈线性增加, 更适用于房屋建筑结构的倒塌仿真分析。閣东东等^[2]在ABAQUS平台上进行二次开发, 采用显式算法对框架结构、框架-剪力墙结构进行了倒塌仿真模拟。

　　 本文基于ABAQUS软件显式计算平台进行二次开发, 自行编制适用于倒塌数值仿真的计算子程序, 对构件有限元模型和失效判别标准进行深入研究, 与文献[2]的主要区别在于:1) 剪力墙的模拟不同, 文献[2]采用带有钢筋层的分层壳, 本文采用弥散杆单元的Rebar Layer分层壳;2) 墙、柱钢筋单元受压破坏模拟不同, 文献[2]假定混凝土压坏之后钢筋仍可以继续受压, 直到压应变为0.01才发生压曲失效, 本文设定混凝土压坏退出工作后, 钢筋作为细长杆同时压曲, 退出工作。本文侧重于论述基于倒塌分析的工程应用, 详细的子程序开发原理和过程另行撰文阐述。

　　 通过逐步加大地震波的峰值加速度, 对建筑结构进行不同地震水平作用下的全过程倒塌数值仿真研究。本文首先对某工程案例计算结果可信性进行验证。其次, 准确发现板柱-剪力墙结构抗震薄弱部位, 研究其倒塌机制, 提出具体加强措施, 通过对比计算验证提高抗震性能方法的有效性。

　　 基于ABAQUS显式计算平台自行研发适用于建筑结构在地震作用下的倒塌数值仿真子程序, 采用通用接触算法分析构件之间的接触、碰撞和构件碎片堆积的影响。构件接触面的摩擦力由库伦摩擦分析, 设置罚函数考虑构件之间可能发生的微小滑移量^[3]。自编子程序实现“单元生死”技术, 真实模拟构件的破坏。

　　 采用分层壳 (S4R, S3R单元) 模拟楼板和剪力墙, 分层壳内根据实际设计的配筋设置钢筋网以精确模拟剪力墙的墙身分布筋和楼板钢筋。采用可以考虑杆件剪切变形的基于Timoshenko理论的B31单元模拟梁、柱单元。壳元模型选用ABAQUS自带的混凝土损伤塑性模型 (Concrete Damaged Plasticity Model) , 壳元中分布筋及梁、柱的材料采用自编的Vumat材料, 其混凝土材料应力-应变骨架曲线及加卸载规则参考文献[4]附录C (图1) , 钢筋材料采用随动强化模型 (图2) 。

　　 钢筋压曲失稳应变对结构倒塌概率的影响较大^[5], 地震作用下, 剪力墙通常发生剪压破坏, 故对于含剪力墙构件的结构, 其钢筋受压失稳应变越大, 倒塌概率越小。钢筋混凝土构件中混凝土压碎之后, 钢筋成为失去侧向支撑的细长压杆, 其失稳应变由约束条件、初始缺陷 (侧移) 决定。强震下构件发生较大侧移, 偏于保守考虑, 本文剪力墙、框架柱钢筋受压失稳应变取0.003 3 (混凝土极限压应变) , 墙、柱失效以混凝土压坏或者钢筋拉断为判别标准。考虑到梁和楼板单元受到其周围板单元的约束, 梁和楼板的混凝土压碎之后, 在大变形下其钢筋仍能起到对相连的竖向构件的拉结作用, 因此楼板和梁失效以钢筋拉断^[1]作为失效判别准则。构件失效判别标准详见表1。

　　 图1 混凝土单轴本构模型 (C30)

　　 为验证所研发倒塌计算程序的合理性, 以某板柱-剪力墙结构为例进行研究分析。

　　 项目位于深圳市福田区, 建筑结构安全等级为二级, 设计使用年限为50年, 抗震设防分类为丙类。抗震设防烈度为7度 (0.1g) , Ⅱ类场地。地下3层, 地上13层, 结构总高度为49.8m, 标准层层高为3.7m。结构体系最终采用板柱-剪力墙结构 (图3) , 项目自身特点归纳如下:1) 塔楼总高受规划要点限制, 楼盖体系采用无梁楼盖, 可以增加一个楼层, 提升项目的容积率;2) 楼盖体系与宽扁梁方案比较, 楼盖高度由450mm减小到260mm, 增加楼层净高190mm, 满足建筑对楼层净高要求。

　　 基于建筑平面布置和功能需求, 塔楼2层大部分楼盖缺失 (图4) , 核心筒区域墙体和框架柱形成穿层墙和穿层柱。根据文献[6-7]有关板柱-剪力墙结构的构造要求, 标准层楼盖体系边跨设置框架梁 (图5) 。无梁楼盖楼板厚均为260mm, 框架柱之间按照构造要求设置暗梁, 柱帽高度为500~650mm (图6) 。

　　 根据《超限高层建筑工程抗震设防专项审查技术要点》 (建质[2015]67号) ^[8], 本项目可不计入超限高层建筑范畴, 但是考虑到此类型结构体系受力机理复杂, 国内工程设计经验和震害案例较少, 认为有必要进行性能化设计, 为揭示其抗震性能机理做更为深入的基础性研究。

　　 分别从构件和整体层面验证开发的倒塌计算子程序计算结果的准确性。下文中所述中震、大震、超大震弹塑性时程分析以及倒塌分析均采用自行研发的倒塌计算子程序。采用El Centro波进行结构动力弹塑性时程分析。地震波按文献[5]的要求将峰值加速度调整为220 (7度大震) , 400 (8度大震) , 620cm/s² (9度大震) , 地震波的加速度时程曲线见图7。

　　 以9度大震弹塑性时程分析为例进行构件层面计算结果的可信性验证。图8为底层核心筒剪力墙Wa (图4) 计算结果反演的混凝土等效塑性应变-受压损伤曲线图, 当混凝土塑性应变达到0.02时, 混凝土受压损伤系数为0.977, 混凝土材料失效, 构件发生破坏。图9为剪力墙Wa的混凝土等效塑性应变时程曲线图, 当t=4.2s, 由于其相连单元因失效退出工作, 结构内力重分配, 导致剪力墙Wa等效塑性应变随之突然增大至0.02, 混凝土材料失效, 构件丧失承载力, 与表1构件失效判别准则吻合。

　　 图10为底层框架柱KZa (图4) 计算结果反演的混凝土滞回曲线和文献[4]混凝土应力-应变关系对比图, 混凝土压应变为0.003 3时, 材料失效, 构件发生受压破坏。图11为框架柱KZa计算结果反演的钢筋应力-应变滞回曲线, 结合图10, 11可知, 混凝土压应变为0.003 3时, 材料失效, 钢筋随之发生压曲, 失去承载力, 与表1构件失效判别准则吻合。典型构件的反演表明, 所编写程序计算可靠准确, 能真实仿真模拟钢筋混凝土构件从加载直至破坏的全受力过程。

　　 以7度大震弹塑性时程分析为例, 并同步采用大震弹性时程方法进行相同地震波 (El Centro波) 作用下的对比分析。两种计算方法所提取的基底剪力时程曲线 (图12) 发展规律归纳如下:1) 在0~3s, 大震弹性时程分析与大震弹塑性时程分析的基底剪力基本一致, 曲线基本吻合;2) 在3~6s, 大震弹塑性时程分析的基底剪力逐渐小于大震弹性时程分析的基底剪力, 这是由于随着地震剪力增大, 部分耗能构件 (如:连梁、框架梁、柱帽区域) 屈服, 结构刚度退化;3) 6s以后, 两者基底剪力时程曲线差异逐渐加大, 大震弹塑性时程计算的基底剪力明显小于大震弹性计算的基底剪力, 且两者相位差明显。大震弹塑性时程分析基底剪力最大值和大震弹性时程分析基底剪力最大值的比值为0.54。

　　 相同地震波作用下大震弹塑性时程分析与小震弹性时程分析的基底剪力之比为3.42 (表2) 。按结构在多遇地震、设防地震和罕遇地震的最大地震影响系数比值评估, 该比值位于2.875~6.25区间, 计算结果反映结构进入了弹塑性受力状态, 大震弹塑性计算的基底剪力介于中震弹性和大震弹性之间, 符合抗震受力概念。

　　 图13为El Centro波大震弹塑性时程分析和大震弹性时程分析的顶点位移对比曲线, 由图可见:1) 在0~3s, 大震弹塑性时程分析和大震弹性时程分析的顶点位移基本一致, 曲线基本吻合;2) 3s以后, 耗能构件屈服增多, 结构损伤加大, 周期变长, 大震弹塑性时程分析的顶点位移变小;3) 9s以后, 大震弹塑性时程分析和大震弹性时程分析的顶点位移出现明显的相位差。

　　 由大震弹塑性时程分析和大震弹性时程分析的基底剪力和顶点位移的曲线可见, 在地震波作用下, 两者曲线初始阶段吻合, 到中后期, 由于耗能作用, 存在结构基底剪力、顶点位移逐步衰减的发展过程, 计算结果符合结构抗震概念, 真实可信。

　　 整体和构件层面的验证结果均表明本文开发的倒塌计算程序计算结果符合结构概念, 能反映构件失效的真实情况, 程序可靠性可得到保证。

　　 小震弹性反应谱计算模型的前3阶周期为:T₁=1.95s (Y向平动) , T₂=1.80s (X向平动) , T₃=1.77s (扭转) 。X, Y向剪重比分别为2.07%, 1.95%。周期比和剪重比均满足规范要求。楼层侧向刚度比均大于1 (图14 (a) ) , 属于侧向规则。X向楼层扭转位移比小于1.2, Y向楼层扭转位移比最大值为1.32, 小于1.4 (图14 (b) ) , 属于扭转不规则。由各项指标可见, 结构宏观抗震性能均在可控范围之内。

　　 采用YJK软件对结构进行中震作用下的等效弹性分析, 连梁刚度折减系数取0.4、周期折减系数取1.0、阻尼比取5%。参考文献[7], 中震作用下竖向构件性能目标为抗剪弹性、抗弯不屈服。计算结果显示, 结构在底层有较多X向墙肢和部分Y向小墙肢受力状态为偏心受拉, 拉力和墙肢混凝土抗拉承载力标准值之比 (N_t/f_tkA_c) 介于1.5~4.3之间, 见图15, 图中填充构件为偏心受拉墙肢。

　　 基于中震不屈服结果计算剪力墙暗柱配筋, 偏心受拉墙肢的配筋率较大, 暗柱最大配筋率为6.17% (图16) 。按中震不屈服计算结果进行试设计, 以暗柱AZ1为例, 其纵筋配筋率为5.25%, 配筋详图见图17, 纵筋平均净间距仅为45mm, 暗柱纵筋密集, 其设计结果是否合理和经济有待商榷。

　　 文献[9]第4条第3款指出:小偏拉钢筋混凝土构件混凝土开裂、钢筋屈服后, 其侧向刚度、轴向受拉刚度明显下降, 所分配的拉力、水平剪力明显减小。为简化设计, 可偏安全地按等效弹性计算分配的拉力、水平剪力进行构件截面承载力计算, 并宜按中震弹性或不屈服设计, 必要时也可采用弹塑性方法进行内力分析。基于上述成果, 采用ABAQUS软件建立弹塑性模型, 用中震水平的El Centro波 (PGA=100cm/s²) 加载计算。分别提取YJK软件中震等效弹性方法计算所得的墙肢W1 (图15) 拉力和ABAQUS中震弹塑性时程分析提取的墙肢轴力时程 (图18) 。前者墙肢W1拉力值为3 237k N, 后者墙肢W1最大拉力值为2 390k N。两者墙肢拉力之比为0.74。上述对比分析表明中震不屈服算法偏于保守, 墙肢计算内力偏大, 中震弹塑性时程分析结果更接近墙肢的真实受力状态。

　　 塔楼设计先后经历两轮, 原结构经过初步计算后发现抗震能力较弱, 经过采取加强措施后予以重新复核, 抗震性能得到明显改善, 下文进行详细的对比计算分析。

　　 对原结构进行在设防大震 (7度, PGA=220cm/s²) 作用下的弹塑性时程分析。图19为剪力墙损伤云图, 墙肢损伤总结如下:1) X向连梁两端的墙肢沿塔楼高度损伤分布均匀, 受压损伤系数介于0.1~0.3之间;2) 底层X向部分墙肢有较严重损伤, 受压损伤系数为0.87;3) 个别墙肢在底部发生轻微弯曲屈服, 其他部位损伤较小, 沿塔楼高度均匀分布, 起到耗能作用。

　　 图1 9 原结构7度大震作用下剪力墙损伤云图

　　 图20为框架梁、柱损坏等级分布图, 塔楼中上部框架柱损坏等级为4, 为中度损坏。连梁、框架梁损伤等级分别为3和2, 为轻度和轻微损坏, 损坏部位沿高度均匀分布, 起到耗能作用。

　　 图2 0 原结构7度大震作用下框架梁、柱损坏等级分布图

　　 图21为8层楼板受压损伤云图, 楼板在核心筒连梁和框架柱的柱帽位置损伤较大, 损伤系数达到0.52, 其余位置楼板损伤较小, 处于基本完好状态。

　　 图2 1 原结构7度大震作用下8层楼板损伤云图

　　 图22为8层楼板钢筋最大受拉应变云图, 应变超过0.002的区域表示楼板内钢筋出现屈服。楼板钢筋在核心筒连梁位置及框架柱柱帽的顶部区域发生轻度屈服, 最大拉应变为0.003 6, 其余位置楼板钢筋保持弹性。柱帽区域出现屈服, 表明板柱-剪力墙结构由于自身体系的缺陷导致其抗震性能相对较弱。

　　 图2 2 原结构7度大震作用下8层楼板钢筋拉应变云图

　　 综上所述, 原结构在设防大震下未发生倒塌, 但是底层部分墙肢损伤严重, 无梁楼盖柱帽区域普遍发生轻微屈服。

　　 原结构在加载8度大震t=6.5s时发生倒塌 (图23) , 关键时刻如下:1) t=4s时, 连梁、柱帽外边缘出铰, 底部加强区墙体、框架柱均未出铰;2) t=4.5s时, 框架梁、柱帽大量出铰, 结构右侧墙体剪力墙底部加强区暗柱钢筋拉断, 开始出现剪压破坏, 混凝土压碎退出工作;3) t=5s时, 框架梁出铰进一步增加, 结构右侧筒体剪力墙破坏进一步加剧, 基底剪力向左侧筒体转移, 左侧筒体底部加强区破坏严重, 同时右侧部分框架柱出铰并退出工作, 结构出现倒塌趋势;4) t=5.5s时, 底部两层框架柱全部出铰并退出工作, 底部两层结构倒塌, 引起整体结构向右侧倾斜;5) t=6.5s时, 由于上部结构的撞击, 整体结构垮塌严重。

　　 基于上述初步研究成果, 为提高结构抗倒塌能力, 经多轮反复计算, 提出如下改进意见:1) 底部1~2层核心筒外周圈剪力墙 (图24中的填充墙体) 墙厚加厚至400mm;2) 上述加厚剪力墙暗柱配筋率提高至2%, 墙身配筋率提高至0.4%, 其他不变。

　　 若暗柱配筋率取2%, 经试设计分析, 暗柱钢筋间距合理, 其实际配筋符合结构常规概念。以AZ1 (图16) 为例, 其详图见图25。

　　 原结构和改进后结构的ABAQUS大震弹塑性层间位移角曲线如图26所示, 在7度大震作用下, 最大层间位移角由1/236减小为1/275。在8度大震作用下, 改进后结构最大层间位移角为1/151, 结构未发生倒塌。通过针对性优化, 塔楼抗侧刚度有所提高, 即便是在8度超大震作用下也可满足性能目标要求。通过精细化的计算分析, 取小于常规中震等效弹性模型计算的剪力墙暗柱配筋值进行设计, 仍可满足“大震不倒”性能目标。

　　 根据文献[6-7]中有关暗梁构造措施的要求, 在核心筒与外框梁之间的无梁楼盖柱上板带区域设置暗梁, 暗梁高度取楼板厚度, 宽度取柱宽及柱两侧各1.5倍板厚之和 (图5) , 暗梁支座上部钢筋截面面积取不小于柱上板带钢筋截面面积的50%, 并全跨拉通, 暗梁下部钢筋截面面积取不小于上部钢筋截面面积的50%。

　　 为考察无梁楼盖中暗梁的作用, 建立不设置暗梁的模型并进行8度大震弹塑性时程分析, 计算得到最大层间位移角为1/140 (图26) , 地震波加载结束时, 结构未发生倒塌。与设置暗梁的模型相比, 1~7层的最大层间位移角有所增加, 最大层间位移角增加了8%。表明暗梁参与了耗能, 减小了结构的损伤程度, 对提高结构抗侧刚度和抗倒塌能力有一定帮助。

　　 8度大震作用下, 不设置暗梁模型与设置暗梁模型楼板损伤区域差异较大 (图27) , 具体表现为:1) 在不设置暗梁模型中, 核心筒两侧的大跨楼板区域中度损伤 (损伤系数为0.63) , 设置暗梁后该区域变为基本完好;2) 在不设置暗梁模型中Y向柱跨的柱上板带和跨中板带损伤严重 (受压损伤系数0.87) , 设置暗梁后这两个区域表现为轻微损坏~轻度损坏。上述楼板在重力荷载作用下均表现为基本完好, 所发生损伤为地震往复作用造成, 表明暗梁在地震往复作用下可有效耗散能量, 进而起到减小楼盖结构损伤的作用。同时, 上述结果验证了规范^[6]对设置暗梁相关规定的必要性。

　　 从顶点位移曲线和构件两个层面分别研究结构的抗震性能。由顶点位移曲线可发现结构整体损伤规律, 构件损伤的研究能从微观上反映结构抗震薄弱部位。

　　 图28为原结构在7, 8度大震作用下弹性、弹塑性时程分析顶点位移时程曲线。7度大震作用下, 弹性时程分析和弹塑性时程分析的结构顶点位移以水平轴为中心往复摆动, 两者曲线基本同步、吻合, 表明构件对变形具有良好的恢复特性。8度大震作用下, t=5s前, 顶点位移以水平轴为中心往复摆动, t=5s后, 结构顶点位移值呈直线单调增加, 表明结构发生不可逆转的侧向倾斜, 直至倒塌。

　　 图29为加强后结构在7~9度大震作用下的弹塑性时程分析顶点位移时程曲线。由图可见:1) 7, 8度大震作用下, 结构顶点位移以水平轴为中心往复摆动, 表明构件对变形具有良好的恢复特性;2) 9度大震作用下, t=5s前, 顶点位移以水平轴为中心往复摆动, t=5s后, 结构顶点位移值呈直线单调增加, 表明结构发生不可逆转的侧向倾斜, 直至倒塌;3) 随着地震烈度的增加, 顶点位移曲线摆动幅度增大, 相位差逐步加剧, 表明结构损伤构件增多, 且自振周期变长;4) 对比图28, 29中顶点位移曲线可见, 采取适当的加强措施后结构体系抗倒塌能力相当于提高1度。

　　 加强后结构在9度大震 (El Centro波) 作用下的倒塌计算结果如图30~33所示。由图可见:1) t=4s时, 核心筒连梁、外框梁、柱帽外边缘沿塔楼高度均匀屈服, 结构底部个别剪力墙发生严重损伤, 中上部墙体发生中度损伤, 框架柱发生中度损坏 (图30) ;2) t=4.5s时, 梁及柱帽大量出铰, 右侧筒体底层的墙肢W11 (图24) 发生剪压破坏, 随后靠右侧墙肢W23, W12, W25发生剪压破坏 (图31) ;3) t=5s时, 墙肢W2, W3, W4, W8, W9, W10发生剪压破坏, 墙肢破坏逐步向左侧发展, 右侧部分框架柱出铰并退出工作, 结构发生倾斜, 出现倒塌趋势 (图32) ;4) t=5.5s时, 底部剪力墙均发生破坏, 框架柱基本出铰并退出工作, 结构底层发生倒塌, 引起结构整体向右侧倾斜倒塌 (图33) 。

　　 图34为另一组9度大震水平地震波 (人工波) 作用下结构的倒塌计算结果。由图可知:1) t=4.5s时, 结构部分连梁、框架梁、柱帽外边缘出铰;2) t=5s时, 梁铰、柱帽铰进一步增加, 结构右侧的筒体底部加强区出现塑性铰;3) t=5.5s时, 右侧筒体底层剪力墙破坏, 右侧部分框架柱出铰;4) t=6s时, 左侧筒体底层剪力墙剪压破坏, 结构整体呈向右垮塌趋势, 框架梁、楼板的拉结作用明显;5) t=7s时, 底层结构垮塌。

　　 由以上两组地震波的倒塌分析过程可见, 本结构墙柱的塑性铰最先出现在底部加强区, 并且由于该底层墙柱的混凝土压碎从而退出工作, 导致整体结构最终倒塌, 两条地震波倒塌破坏次序在宏观层面上基本一致, 更进一步验证了计算结果的可靠性。

　　 基于实际工程的需求, 在ABAQUS软件显式计算平台上, 开发适用于模拟地震作用下建筑结构倒塌的计算程序, 并从整体和构件层面验证程序计算结果准确性。通过中震、大震、超大震计算以及倒塌分析, 研究板柱-剪力墙结构抗震性能, 研究成果汇总如下:

　　 (1) 有别于梁板结构, 无梁楼盖体系柱帽区域在大震作用下出现轻度屈服, 有必要严格执行文献[6-7]中有关板柱-剪力墙结构设置暗梁的构造措施。

　　 (2) 通过精细化的中震弹塑性分析验证, 等效弹性算法计算的中震不屈服结果较为粗略, 中震弹塑性分析更接近实际受力状态, 且经济性优势明显。

　　 (3) 全过程倒塌分析表明, 适当提高底部加强区墙体和暗柱配筋率对结构抗倒塌能力有显著改善。

　　 (4) 针对板柱-剪力墙结构, 在常规设计方法存在不合理情况时, 为兼顾结构体系安全性与经济性, 有必要进行中震、大震弹塑性时程和倒塌分析, 做到有针对性的加强结构设计。