随着城市规模的扩大、人口密度的增加以及现代工业的规模化发展，地震作用下建筑物严重破坏和倒塌导致的人员伤亡和经济损失越来越大，而钢筋混凝土框架结构是目前我国建筑中应用最为广泛的一种结构形式，因此框架结构的抗震性能和抗倒塌性能研究也就显得格外迫切。
关于钢筋混凝土框架结构地震下的倒塌研究，各国学者已经取得了一定的成果。William[1]比较了不同单元的优缺点：塑性铰杆单元模型计算最为简单快速，适合用于大型结构计算；纤维单元模拟精度最高；三维实体模型最为复杂，计算往往不易收敛。曹鹏等[2]基于能量等效假设开发了用于ABAQUS实体单元的混凝土塑性损伤模型。周浩等[3]采用三种混凝土本构关系模拟结构倒塌，对比说明选取适当本构关系的重要性。冯德成等[4]基于显式理论提出了KR-α法模拟结构倒塌，兼顾计算效率和精度的同时，时间步长可以更大。郭金龙等[5]模拟了地震作用下框架结构的落层碰撞效应。陆新征课题组[6－8]基于MARC平台，开发了能模拟钢筋混凝土杆系结构的纤维模型程序THUFIBER，对地震作用下框架结构倒塌数值仿真进行了系统的研究，研究表明高层结构随着输入地震动的不同，结构的倒塌发生楼层和倒塌模式会出现很大变化，但薄弱层主要位于下部楼层。顾祥林等[9]基于离散单元法，开发了钢筋混凝土框架结构空间倒塌反应数值仿真系统，并通过相应的振动台试验验证了该系统的正确性。强震作用下结构的倒塌破坏属于大变形不连续问题，传统的有限元隐式积分通常难以满足其收敛条件，离散元法虽能满足位移不连续性，但分析精度较低，耗时长[10]。
本文在以上研究的基础上，基于 ABAQUS显式分析平台，开发了可以考虑应变率效应及强度和刚度退化效应的三维梁单元材料子程序，基于构件损伤定义层损伤，对某六层钢筋混凝土框架结构进行增量动力时程分析，研究其倒塌机理，模拟倒塌过程。
1 基于ABAQUS平台的本构开发
1.1 关于ABAQUS隐式和显式模块
对于ABAQUS 平台，有隐式和显式两大求解模块。隐式方法在静力分析和弹性动力时程分析中，求解快，精度高；但在高度非线性和大型结构计算中，需要不断缩小步长迭代求解，每次迭代都要重组刚度矩阵，求解代价大，且极容易计算不收敛。显式方法采用中心差分法，应用前一个增量步的动力学条件计算下一个增量步的动力学条件，它的求解依赖于一个稳定时间步长，不存在迭代收敛问题，对于某些大型非线性问题会更具有优势。
倒塌分析中结构单元数目较多，且存在大变形和高度非线性的问题，如用隐式算法效率较低，还可能不收敛，因此选用显式算法比较合适。
1.2 钢筋材料本构模型开发
虽然构件承载力退化是构件的宏观行为，但其退化机理通常认为是反复荷载作用下钢筋与混凝土间的粘结滑移和混凝土保护层的剥落导致。同时，钢筋的本构关系对钢筋混凝土结构的杆系构件在往复加载下的压弯行为有决定性的影响。已有学者在研究构件承载力退化行为时，将构件承载力的退化全部计入钢筋纤维的滞回本构中[11]，其与试验结果的比较表明，这种处理方法能有效地模拟钢筋混凝土构件在逐级递增加载下的承载力退化。这说明，将宏观构件的承载力退化计入钢筋的滞回本构是行之有效的。
本文对文献[12]提出的基于有效累积滞回耗能的钢筋混凝土构件承载力退化模型进行优化，开发钢筋的单轴本构关系。该模型在钢筋与混凝土分离模型的基础上，在钢筋恢复力模型中考虑了粘结滑移和混凝土压溃引起的构件承载力退化，钢筋屈服前认为是理想弹性，钢筋屈服后考虑屈服强度的退化，再加载时指向历史最大应变点，且该指向点对应的钢筋屈服强度随累积耗能的增大而减小。这里的强度退化并非指钢筋本身的强度退化，而是反映了钢筋与混凝土之间由于粘结滑移以及混凝土保护层剥落所产生的综合退化效果，如图1所示。图中的σy为初始屈服强度；E0为初始弹性模量，代表钢筋的弹性加载或卸载下降段；εi、Eqi分别为第i次卸载时的应变和反向再加载时相应的等效弹性模量。
图1 钢筋应力-应变曲线
Fig.1 Stress-strain curve of rebar
其中有效累积滞回耗能定义为：
式中：Eeff为当前的有效累积滞回耗能；ΔE表示当前增量应变下的有效能量；ε为当前应变；εy为钢筋的屈服应变；εmax为历史最大应变。
钢筋骨架线加载关系如下式所示，为理想弹塑性骨架：
式中：σ为当前应力；Es为初始弹性模量，根据规范按照钢筋标号选取；εu为钢筋极限应变，可以按照相关试验数据选取，也可以按照相关经验公式计算；σy为初始钢筋屈服强度，按规范取标准值或设计值。
钢筋卸载关系如下式所示，为弹性：
式中，σun、εun为卸载点的应力、应变。钢筋反向再加载关系如下式所示：
式中：σyi为折减后的当前屈服强度；Eeff为按照式(1)求得的有效累积滞回耗能；σy为初始钢筋屈服强度；εf可以理解为钢筋与混凝土之间粘结滑移失效时的应变，但其影响因素较为复杂，为简便起见，根据模拟结果与试验结果对比，取εf=25εy时拟合情况较好，其中εy为钢筋的屈服应变；λ、c为通过试验确定的参数，这里根据试验数据标定为λ=6，c=0.3。
式中：Eq为反向再加载时的等效弹性模量；σlast、εlast为上一时刻的应力、应变。
1.3 混凝土材料本构模型开发
本文采用《混凝土结构设计规范》中提出的混凝土本构模型[13]，如图 2所示(受拉及受压曲线均绘于图中，但所取比例不同)。卸载及再加载时均视为弹性，受压状态下从卸载点卸载至残余应变εz处，受拉状态下从卸载点卸载至初始受拉点处，并设置极限压应变εcu和开裂应变εf。
图2 混凝土应力-应变曲线
Fig.2 Stress-strain curve of concrete
1.4 应变率效应
地震作用下，钢筋的应变率最大可达1/s，混凝土的应变率最大可达 10-2s-1[14]，随着应变率的提高，钢筋屈服强度和混凝土的抗压强度有明显的增加，但应变率对于材料的损伤阈值及其与临界应变的比值影响不大，从而导致退化效应更加明显。考虑到计算效率，这里不考虑混凝土弹性模量随应变率的变化。
本文采用大连理工大学李敏提出的钢筋动态本构模型[15]，其屈服强度动态增大系数cf按照下式计算：
式中：σyd和σys分别为钢筋动态和静态屈服强度；
为当前应变率；
为准静态应变率，这里取为
采用大连理工大学肖诗云提出的混凝土动态本构模型[16]，其抗压强度动态增大系数按照下式计算：
式中：fcd和fcs分别为混凝土动态和静态抗压强度，
为当前应变率；
为准静态应变率，这里取为
1.5 结构层损伤定义
结构在地震下发生倒塌破坏的过程其实是损伤发生和不断累积的过程，本文从材料本构出发，将宏观构件的损伤退化计入钢筋的滞回关系，并以钢筋单元损伤来表征结构的层损伤，定义每个单元的损伤值为：
式中：
为第j层第i个构件中单元的损伤指数；E、E0为材料当前等效弹性模量、初始弹性模量；
为该构件的损伤指数，取为构件中单元损伤最大值。
考虑到损伤程度越严重的构件对整体结构的层损伤贡献越大[17]，将层损伤定义为：
式中：Dj为第j层损伤值；n为第j层构件数；mij为损伤加权系数；Dij为式(11)中定义的第j层第i个构件的损伤指数。
1.6 子程序试验验证
本文对文献[18]中的钢筋混凝土柱在静力加载和动力加载两种工况下的受力情况进行模拟，试件构造、配筋及加载方式如图3所示，图中所示尺寸量纲均为mm，纵向钢筋为HRB335，实测屈服强度 381.5 MPa，混凝土为C30，实测抗压强度26.4 MPa。加载规则如图4所示，沿构件的轴线方向施加恒定的轴力，控制轴压比为0.1，侧向施加位移控制的低周往复荷载，第一周加载幅值为5 mm且仅循环一周，接下来的位移幅值以5 mm作为增量直至破坏，每级加载幅值循环三周，整个加载过程中，保持速率恒定，静力加载时，加载速率为0.1 mm/s，动力加载时，加载速率为40 mm/s。
图3 试件构造、配筋及加载方式
Fig.3 Specimen construction, reinforcement and loading
图4 加载规则
Fig.4 Loading rules
图5 静力试验结果与数值模拟对比
Fig.5 Comparison of experimental results and numerical simulations under static
图6 动力试验结果与数值模拟对比
Fig.6 Comparison of experimental results and numerical simulations under dynamic
试验与模拟结果对比如图5和图6，可以看出无论对静力试验还是动力试验，在 ABAQUS平台中使用本文所开发的材料本构模型，均能比较准确的描述构件的滞回行为，特别是很好地模拟了构件的承载力退化行为。动静态的滞回曲线也体现出材料的应变率效应：构件的峰值承载力有一定提高，但同时退化现象也更加明显。因此，地震作用下的结构分析有必要考虑应变率效应和累积损伤退化效应。
2 框架结构设计与建模
结构选用某六层钢筋混凝土框架结构，使用PKPM (2010规范)计算配筋，结构的模型尺寸如图7所示，框架层高为3.6 m，底层层高考虑基础的埋深为950 mm，框架梁在x方向的截面尺寸取为300 mm×600 mm，在y方向的截面尺寸取为300 mm×700 mm，次梁的截面尺寸取为250 mm×500 mm，框架柱的截面尺寸均取为600 mm×600 mm，相关的设计信息如表1所示，材料选用标号HRB335的钢筋及C30混凝土，具体属性见表2、表3。
图7 结构示意图 /mm
Fig.7 Schematic diagram of structure
表1 设计信息
Table 1 Structure design information
表2 钢筋属性
Table 2 Mechanical properties of HRB335
表3 混凝土属性
Table 3 Mechanical properties of C30
使用 ABAQUS模拟时考虑重力荷载代表值，将构件自重以外的荷载均匀地添加到楼板上。梁柱构件均使用*ELCOPY命令进行配筋，采用B31单元，材料采用上面编制的子程序。研究表明在倒塌分析中楼板对梁的约束应予以考虑[19]，这里取楼板厚度为100mm，顶层120mm，采用S4R单元，设置为弹性，不考虑楼板的损伤和率效应。有限元模型如图8所示。
3 结构倒塌模拟及分析
图8 结构三维有限元模型
Fig.8 Three-dimensional finite element model of structure
对结构进行模态分析，其一阶振型为x方向平动，周期为1.2617 s，因此选择x方向为地震施加主方向，并按照规范，在其他方向按照 1∶0.85∶0.65的比例进行加载。
模拟时采用生死单元法，当材料达到破坏准则时，移除单元，模拟破坏，这里的破坏准则设置为：当钢筋应变|ε|超过极限应变εu时，钢筋视为断裂，单元不再具有承载力，从模型中移除；混凝土的破坏准则取为当压应变ε＞εcu时，混凝土视为压碎，单元被移除，对于受拉情况，则设置开裂应变，当混凝土拉应变超过开裂应变后，混凝土应力恒定为0，不再具有抗拉能力，但并不移除单元，仍具有受压承载力。
选用EL Centro波，对结构进行增量动力时程分析(IDA)，将地震波峰值加速度调幅到1 m/s2，加速度增量设置为0.5 m/s2，对结构进行分析，得到相关结构响应数据，直至结构发生倒塌破坏。
3.1 结构响应分析
图9和图10分别给出了在不同峰值加速度地震作用下结构顶层位移时程曲线和首层层间位移角时程曲线，图 11给出了结构的能力曲线。可以看出，x和y向的位移峰值和层间位移角峰值随地震强度的增大而增大，当峰值加速度达到0.3g时，结构的位移时程曲线和层间位移角时程曲线随时间变化趋于无穷大，结构在地震持时内发生了倒塌破坏。
图9 顶点位移时程曲线
Fig.9 Top displacement time history curve
图10 底层位移角时程曲线
Fig.10 Displacement angle time history curve of first floor
图11 结构能力曲线
Fig.11 Capacity curve of the structure
图 12给出了在不同峰值加速度地震作用下层间最大位移随楼层的变化曲线，可以看出，层间位移角在1层～2层随地震强度增大而急剧增加，其他各层层间最大位移角虽有增加，但幅度基本一致，也就是说该框架结构底部两层是薄弱层，在极限荷载作用下最先发生损伤破坏。
图12 结构最大层间位移角
Fig.12 Maximum inter story drift
图13 顶点位移时程曲线(PGA=0.3g)
Fig.13 Top displacement time history curve (PGA=0.3g)
图14 结构耗能时程曲线(PGA=0.3g)
Fig.14 Structural energy dissipation time history curve(PGA=0.3g)
图13和图14分别给出了在峰值加速度为0.3g的地震作用下，考虑和未考虑损伤累积效应对结构耗能曲线和顶点位移时程曲线的影响。可以看出，当结构处于弹性状态时，耗能和顶点位移曲线基本重合；结构进入塑性状态后，发生损伤并不断累积，在考虑损伤退化的情况下，材料屈服强度下降，应力与应变增量围成的区域相对不考虑损伤退化的情况减小，即耗能能力下降，同时结构变形增大，最终结构发生倒塌破坏后，顶点位移和耗散能量随着时间的变化趋于无穷大，而未考虑损伤退化的结构未发生倒塌破坏。可见，考虑材料损伤和应变率效应后，更准确地表征了结构承载力的退化，这种情况下在地震持时内结构的失效破坏，更符合实际工程结构的退化、倒塌情况。
3.2 倒塌机理分析
图15 结构失效倒塌示意图
Fig.15 Schematic diagram of structural collapse
图16 结构倒塌破坏变形示意图
Fig.16 Structure collapse deformation diagram
图15给出了当峰值加速度为0.3g时所模拟的框架结构倒塌的全过程示意图，图上所示输出变量为单元损伤指数值。图 16为结构倒塌时刻变形示意图。
可以看到，1.5 s以前各构件损伤几乎为0，基本处于弹性阶段；随着地震动的持续，结构开始进入塑性阶段，底部一、二层梁端首先出现相对较大的损伤；随着时间的推移，损伤不断累积，柱中开始出现损伤，8.5 s跨中柱底损伤达到极限，该柱失去承载力，此时结构仍未倒塌，柱上部的梁起到悬索作用，将上部荷载传导至其他柱，9 s时边柱损伤也基本到达极限，底层基本失去承载力，结构即将倒塌。图 17给出了在倒塌情况下，计算得到的框架层损伤发展时程曲线。可以看出，楼层越低其整体损伤程度越大。损伤主要集中在框架底部三层，首层的损伤尤为严重，在地震动峰值加速度出现时段内迅速增加，其他各层出现损伤的时间相对底层有一些滞后，并且其增大的幅度呈逐渐平缓的趋势。9 s时底层损伤达到1，失去承载力，其余各层损伤均迅速提高，二层也基本接近极限。
图17 结构层损伤(PGA=0.3g)
Fig.17 Story damage of structure (PGA=0.3g)
基本上该框架结构的地震倒塌是由于首层框架柱损伤严重，逐渐丧失其竖向承载力，其失效过程为：首层梁端→第二、三层梁端→首层中柱→首层边柱→整体倒塌。所模拟的框架结构发生倒塌破坏的时间为10 s，并未发生在地震加速度峰值处，其倒塌机理在于钢筋混凝土框架结构在地震作用下损伤不断累积，产生了永久损伤，使结构承载能力逐步退化，直至无法承受上部荷载，从而导致整个结构的倒塌。
4 结论
本文基于 ABAQUS平台编制了纤维梁单元材料子程序，可以考虑材料应变率效应以及累积损伤导致的强度刚度退化效应。采用该程序对钢筋混凝土框架结构进行地震作用下的增量动力时程分析，得出以下结论：
(1) 地震作用下结构的应变率较大，累积损伤导致的退化效应显著，在结构分析时应予以考虑。
(2) 考虑材料应变率效应和损伤累积效应能更精确确定框架结构的极限荷载力；采用所编制的程序可追踪结构的损伤，确定结构失效路径和薄弱层。
(3) 框架结构并非在地震峰值强度处发生倒塌，而是由于损伤不断累积，结构逐渐失去承载力而倒塌。
参考文献：
[1]Cofer W F.Documentation of strengths and weaknesses of current computer analysis methods for seismic performance of reinforced concrete members [R].University of California, 1998: 37―38.
[2]曹鹏, 冯德成, 沈新普, 等.基于 ABAQUS平台的塑性损伤子程序开发及其稳定性研究[J].工程力学,2012, 29(SupplⅡ): 101―106.Cao Peng, Feng Decheng, Shen Xinpu, et al.Development of plasticity-damage model based on ABAQUS and algorithm stability analysis [J].Engineering Mechanics, 2012, 29(SupplⅡ): 101―106.(in Chinese)
[3]周浩, 李杰.基于不同本构模型的混凝土结构地震倒塌对比分析[J].建筑结构学报, 2016, 37(9): 8―18.Zhou Hao, Li Jie.Comparison study on seismic collapse of concrete structures with different constitutive models[J].Journal of Building Structures, 2016, 37(9): 8―18.(in Chinese)
[4]Feng D, Kolay C, Ricles J M, et al.Collapse simulation of reinforced concrete frame structures [J].Structural Design of Tall & Special Buildings, 2016, 25(12): 578―601.
[5]郭金龙, 蔡健, 陈庆军, 等.地震作用下钢框架落层倒塌效应模拟分析[J].工程力学, 2016, 33(6): 171―179.Guo Jinlong, Cai Jian, Chen Qingjun.Simulation analysis on pancake collapse response of a steel frame under seismic action [J].Engineering Mechanics, 2016,33(6): 171―179.(in Chinese)
[6]陆新征, 林旭川, 叶列平, 等.地震下高层建筑连续倒塌数值模型研究[J].工程力学, 2010, 27(11): 64―70.Lu Xinzheng, Lin Xuchuan, Ye Lieping.Numerical models for earthquake induced progressive collapse of high-rise buildings [J].Engineering Mechanics, 2010,27(11): 64―70.(in Chinese)
[7]Lu Xiao, Lu Xinzheng, Guan Hong, et al.Collapse simulation of reinforced concrete high-rise building induced by extreme earthquakes [J].Earthquake Engineering & Structural Dynamics, 2013, 42(5): 705―723.
[8]Li Y, Lu X, Guan H, et al.An energy-based assessment on dynamic amplification factor for linear static analysis in progressive collapse design of ductile RC frame structures [J].Advances in Structural Engineering, 2015,17(8): 1217―1225.
[9]顾祥林, 黄庆华, 汪小林, 等.地震中钢筋混凝土框架结构倒塌反应的试验研究与数值仿真[J].土木工程学报, 2012, 45(9): 36―45.Gu Xianglin, Huang Qinghua, Wang Xiaolin, et al.Special collapse responses of reinforced concrete frame structures under earthquakes—test study and numerical simulation [J].China Civil Engineering Journal, 2012,45(9): 36―45.(in Chinese)
[10]宣纲, 顾祥林, 吕西林.强震作用下混凝土框架结构倒塌过程的数值分析[J].地震工程与工程振动, 2003,23(6): 24―30.Xuan Gang, Gu Xianglin, Lü Xilin.Numerical analysis of collapse process for RC frame structures subjected to strong earthquake [J].Earthquake Engineering and Engineering Vibration, 2003, 23(6): 24―30.(in Chinese)
[11]Youssef Maged, Ghobarah A.Strength deterioration due to bond slip and concrete crushing in modeling of reinforced concrete members [J].ACI Structural Journal,1999, 96(6): 956―968.
[12]曲哲, 叶列平.基于有效累积滞回耗能的钢筋混凝土构件承载力退化模型[J].工程力学, 2011, 28(6): 45―51.Qu Zhe, Ye Lieping.Strength deterioration model based on effective hysteretic energy dissipation for RC members under cyclic loading [J].Engineering Mechanics, 2011, 28(6): 45―51.(in Chinese)
[13]GB 50010-2010, 混凝土结构设计规范[S].北京: 中国建筑工业出版社, 2010.GB 50010-2010, Code for design of concrete structures[S].Beijing: China Architecture Industry Press, 2010.(in Chinese)
[14]李敏.材料的率相关性对钢筋混凝土结构动力性能的影响[D].大连: 大连理工大学, 2011: 26―55.Li Min.Effects of rate-dependent properties of material on dynamic properties of reinforced concrete structure[D].Dalian: Dalian University of Technology, 2011:26―55.(in Chinese)
[15]李敏, 李宏男.建筑钢筋动态试验及本构模型[J].土木工程学报, 2010, 43(4): 70―75.Li Min, Li Hongnan.Dynamic test and constitutive model for reinforcing steel [J].China Civil Engineering Journal, 2010, 43(4): 70―75.(in Chinese)
[16]肖诗云, 张剑.不同应变率下混凝土受压损伤试验研究[J].土木工程学报, 2010, 43(3): 40―45.Xiao Shiyun, Zhang Jian.Compressive damage experiment of concrete at different strain rates [J].China Civil Engineering Journal, 2010, 43(3): 40―45.(in Chinese)
[17]Park Y J，Ang A H.Mechanistic seismic damage model for reinforced concrete [J].Journal of Structural Engineering, 1985, 111(4): 722―739.
[18]王德斌.考虑动力效应的钢筋混凝土柱抗震性能研究[D].大连: 大连理工大学, 2011: 28―75.Wang Debin.Dynamic effects on seismic behavior of reinforced concrete column[D].Dalian: Dalian University of Technology, 2011: 28―75.(in Chinese)
[19]王少杰, 刘福胜, 徐赵东.RC空间框架结构竖向倒塌全过程试验研究与理论分析[J].工程力学, 2015,32(5): 162―167.Wang Shaojie, Liu Fusheng, Xu Zhaodong.Experimental and theoretical analysis on entire course of vertical progressive collapse of spatial RC frame structure[J].Engineering Mechanics, 2015, 32(11): 162―167.(in Chinese)