赵桂峰， 马玉宏， 陈嘉佳， 黄浩贤

(1.广州大学 土木工程学院，广州 510006； 2.广东省地震工程与应用技术重点实验室，广州 510405；3. 广州大学 工程抗震研究中心，广州 510405)

自从GB 18036—2015《中国地震动参数区划图》颁布以来，极罕遇地震作用的研究得到了普遍的重视。对于消能减震体系，如何在极罕遇地震作用下保证阻尼器性能发挥、确保减震结构的安全性是当前消能减震领域的关键问题。

近年来，为了使阻尼器能够在中、小震时就能够充分发挥作用，阻尼器响应放大技术得到了国内外学者的重视，研发了连杆机构[1-4]、齿轮机构[5]、杠杆机构[6-10]等多种基于不同机械原理的响应放大装置。研究表明：该技术可以使串联的位移相关型或速度相关型阻尼器的位移、速度等响应放大，从而在中、小震作用下就能够充分发挥阻尼器的耗能能力，取得较好的减震效果；同时由于阻尼器的响应放大作用，阻尼器的设置数量可以减少，从而取得显著的经济效益。但是，当遭遇超出预期的极罕遇地震或更大地震作用时，针对罕遇地震作用而设计的现有各类响应放大装置极有可能会使串联阻尼器放大了的位移、阻尼力等响应超过阻尼器的极限能力，使阻尼器性能下降甚至失效。若采用减小响应放大倍数的措施来防止阻尼器失效，则又可能降低中小震作用下的减震效果。因此，现有各类响应放大装置难以同时兼顾中、小震作用下的减震效果与超出预期的极罕遇地震作用下阻尼器不失效两种情况。

针对上述不足，提出了凸轮式响应放大黏滞阻尼器[11](cam-type response amplification viscous damper， CRAD-VD)，如图1所示，并对其进行了伪静力试验与数值模拟分析[12]、单自由度减震体系的控制效果及能量分析[13]，结果表明：该装置通过凸轮式机构的往复运动，使串联阻尼器在取得响应放大效应的同时，有效地避免了阻尼器在各级地震作用下的位移失效问题；通过串联小吨位的阻尼器就可以取得直接安装大吨位阻尼器的减震效果。

图1 CRAD-VD构造示意图Fig.1 Schematic diagram of CRAD-VD

但是，凸轮式响应放大黏滞阻尼器的理论恢复力公式较为复杂，仅依赖所开发MATLAB程序对单自由度或多自由度体系计算分析难以精细地体现三维实体模型的受力和变形，也难以准确模拟构件的损伤程度变化和分布。现有分析软件SAP2000、ABAQUS、PERFORM3D等均无法通过单个单元或者组合单元来体现CRAD-VD恢复力的非线性特点，因此亟需通过二次开发CRAD-VD耗能单元，对相应减震体系三维非线性模型进行分析。

ABAQUS是功能强大的非线性有限元分析软件，提供了大量用户接口子程序来满足解决复杂问题的需求。通过子程序接口，用户可以按照ABAQUS提供的编译格式编写子程序来定义各类非线性单元。陈聪等[14]基于ABAQUS平台，开发了一种剪切型金属阻尼器；葛少平等[15]开发了基于隐式算法的非线性黏弹性阻尼器；李华伟等[16]开发了适用于钢管混凝土中核心混凝土和考虑包辛格效应钢材的单轴滞回本构模型。

本文基于ABAQUS用户单元子程序编译原理，通过软件提供的用户子程序接口，利用FORTRAN语言开发CRAD-VD耗能单元，开发适用于显式求解器的CRAD-VUEL子单元，并通过与MATLAB计算程序进行对比，验证二次开发的正确性。进而，考虑罕遇与极罕遇地震作用，对安装了CRAD-VD的10层装配式混凝土框架结构的减震性能进行有限元分析，验证二次开发程序的正确性和CRAD-VD装置减震控制的有效性。

1 凸轮式响应放大黏滞阻尼器单元

1.1 CRAD-VD恢复力公式

CRAD-VD可通过人字支撑或墙式支撑安装于结构层间。在地震作用下结构产生层间位移时，CRAD-VD通过支撑推动滚珠丝杠发生水平往复运动，带动滚珠螺母和偏心圆盘旋转，利用凸轮机构推动从动矩形框和黏滞阻尼器发生垂直于丝杠运动方向的竖向往复运动，从而耗散地震能量。

CRAD-VD恢复力计算公式为

(1)

本文选用的黏滞阻尼器的力学模型为

(2)

将式(2)代入式(1)，即为凸轮式响应放大黏滞阻尼器的恢复力计算公式。

1.2 基于ABAQUS的CRAD-VD耗能单元设置

(3)

整体坐标系与单元轴向的转换行向量为

T[-Cx,-Cy,-Cz,Cx,Cy,Cz]

(4)

图2 CRAD-VD空间直杆单元Fig.2 CRAD-VD space straight rod unit

2 基于CRAD-VD耗能单元的非线性消能减震体系的求解算法

ABAQUS软件提供了隐式ABAQUS/Standard和显式ABAQUS/Explicit两种求解器，可求解广泛领域的线性和非线性问题。对于安装CRAD-VD的非线性消能减震体系，其滞回曲线存在高度非线性，考虑到地震动的随机性，隐式算法可能出现由于CRAD参数改变、地震波选取不同导致迭代收敛困难的现象，而显式算法则不需迭代和收敛准则，因此本文采用显式算法开展后续研究。

ABAQUS/Explicit显式求解器采用中心差分法对消能减震体系的非线性运动方程在时间域上进行显式积分，基于上一增量步运动平衡方程的解来计算下一增量步的运动状态。参考ABAQUS用户手册[17]提供的中心差分法的蛙跳格式(Leapfrog)，在t与t+Δt时刻交错求解结构的响应。增量步开始时刻t的结构动力方程为

(5)

(6)

(7)

由式(6)可求出增量步开始时刻t的加速度为

(8)

式(8)中质量矩阵为对角阵，不必同时求解联立方程，因此加速度求解比较简单。

ABAQUS/Explicit求解器在计算速度的变化时假定加速度近似为常数，当前增量步中点速度为前一增量步中点速度加上速度变化值

(9)

式中，Δtt，Δtt+Δt分别为前一增量步和当前增量步的步长。则增量步结束时t+Δt的位移为增量步开始时的位移加上速度对时间的积分

(10)

由式(8)增量步开始时计算满足动力学平衡条件的加速度，进而得到后续时间的速度和位移。显然，每一增量步结束时的响应仅依赖该增量步开始时的响应。分析过程中为了使结果收敛，ABAQUS/Explicit自动控制增量步大小，同时为了提高求解效率，求解器总是尽可能选取稳定极限值作为增量步长。

为了保证增量步中加速度近似为常数，进而确保显式算法的计算精度，时间增量必须非常小，因此分析中增量步大大增多，但由于不必同时求解联立方程组，因此每一个增量步的计算成本低。这说明ABAQUS/Explicit显式算法求解复杂非线性动力问题是有条件稳定的，求解精度要求时间步长要小于某个极限值，即最大步长，因此需要设定稳定极限值，否则会造成数值不稳定、求解不收敛的现象。一般来说，稳定极限值取决于结构的最大频率值ωmax，但计算过程中结构的最大频率受各种因素的影响实时变化，难以准确计算。本文将CRAD-VD恢复力子程序中每个单元赋予一个基于稳定性的时间增量步上限，以免分析中出现发散的现象。ABAQUS/Explicit求解器单元-单元估计法来求解稳定增量步长Δtstable的公式为

(11)

式中：Le为单元长度；cd为材料波速。ABAQUS/ Explicit在求解时，若存在形状较小或不规则的单元，为了增大稳定增量步长，提高计算效率，该软件提供质量缩放技术减小材料波速，或放大网格尺寸以增长单元长度，同时保证分析精度，节省分析时间。

依据以上算法，基于VUEL-CRAD耗能单元的非线性消能体系显式算法程序的流程图，如图3所示。主程序调用VUEL子程序，先把初始状态代入运动平衡方程求出初始加速度，再求解下一步的节点位移以及速度，同时利用位移、速度算出地震作用、恢复力、阻尼力等值，进而利用运动平衡方程求解出下一步节点的加速度，以此类推，进行求解。

图3 CRAD-VD消能体系显式算法编写流程图Fig.3 Flowcharts of explicit algorithm for energy dissipation systems with CRAD-VD

3 二次开发程序的可靠性验证

3.1 简谐荷载下二次开发程序的可靠性分析

采用MATLAB软件编写简谐波位移加载的CRAD-VD耗能单元计算程序，并与ABAQUS计算结果进行对比分析。对CRAD-VD耗能单元，约束一个节点的3个平动自由度，对另一节点施加轴向的简谐荷载，加载工况，如表1所示。

表1 简谐波加载工况Tab.1 Simple harmonic displacement loading conditions

ABAQUS显式算法与MATLAB计算的单元滞回曲线对比，如图4所示。可见，滞回曲线吻合很好。在计算过程中，恢复力峰值点处可能会出现微小的抖动，原因是中心差分法基于加速度求解动力响应，而CRAD-VD恢复力的高度非线性特征使其加速度呈非光滑变化。ABAQUS/Explicit显式求解器通过定义体积黏性引入因体积应变造成的阻尼，从而解决了高度非线性可能产生的数值抖动问题[18]。

图4 不同工况下两种软件滞回曲线对比Fig.4 Comparison of hysteretic curves of MATLAB and ABAQUS under different conditions

3.2 地震作用下二次开发程序的可靠性

以CRAD-VD单自由度消能减震体系为例，如图5所示，对比分析地震作用下MATLAB和ABAQUS软件计算的结构响应。单自由度结构的参数为：刚度K=2 000 kN/m，质量M=50 660 kg，阻尼比5%，结构自振周期为1 s。选取1940年El Centro南北向地震动加速度时程，PGA为0.4g。CRAD-VD的设计参数为：阻尼器黏滞系数10 kN·(s/m)0.4，速度指数均为α=0.4，螺距Ld=0.02 m，偏心距e=0.1 m，凸轮半径r=0.1 m，丝杠半径rss=0.025 m，摩擦因数μ=0.1，传动效率η=0.95。

图5 CRAD-VD单自由度减震体系的计算模型Fig.5 Computational model of CRAD-VD single-degree-of-freedom damping system

由MATLAB与ABAQUS显式算法计算的CRAD-VD单自由度减震体系动力反应时程见，如图6所示。CRAD-VD整体滞回曲线及该装置中所串联的黏滞阻尼器C-VD单独的滞回曲线，如图7所示。可见，两种软件的计算结果吻合良好。在此基础上，对不同模型参数和不同输入地震波进行对比分析，均表明所开发的显式算法是可靠的，可直接用于CRAD-VD三维模型的计算分析。

图6 单自由度减震体系反应时程曲线对比Fig.6 Comparison of time curves for SDOF system reaction

图7 CRAD-VD及C-VD滞回曲线Fig.7 The hysteretic curves of CRAD-VD and C-VD

4 装配式消能体系的抗震性能分析

4.1 ABAQUS结构模型建立与地震波选取

以10层装配式混凝土框架结构为研究对象，层高3.9 m，梁柱、楼板混凝土分别采用C35、C30，主筋为HRB400钢筋。抗震设防烈度8度，场地类别Ⅱ类，设计基本地震加速度0.2g，设计地震分组为二组。梁柱采用B31单元，板采用S4R单元。钢筋采用典型的双折线硬化本构模型，混凝土损伤本构采用GB 50010—2010《混凝土结构设计规范》建议的单轴受力本构模型。为了体现装配式结构的特点，对混凝土和钢筋的刚度和承载力进行折减。

在该结构上分别安装凸轮式响应放大黏滞阻尼器CRAD-VD和普通黏滞阻尼器VD，对比研究罕遇与极罕遇地震作用下两种阻尼器对装配式混凝土框架结构的减震控制效果。消能部件宜设置在结构变形较大的位置，并沿两个主轴方向分别设置，形成均匀合理受力体，具体布置如图8所示。

图8 结构平面及阻尼器布置图Fig.8 Layout plan of structure and damper

CRAD-VD采用所开发的VUEL-CRAD显式单元，丝杠螺距Ld=0.03 m，凸轮偏心距e=0.02 m，凸轮半径r=0.1 m，CRAD的设计最大位移、最大速度和最大阻尼力放大倍数分别为2.67，4.18，6.28。设CRAD-VD串联黏滞阻尼器的阻尼系数之和为600 kN·(s/m)0.3，普通黏滞阻尼器VD的阻尼系数为1 500 kN·(s/m)0.3，速度指数均取0.3，二者的阻尼系数之比为0.4。所建立的三维结构空间模型,如图9所示。

图9 基于ABAQUS的分析模型Fig.9 Analysis model based on ABAQUS

采用基于性态的抗震设计思想，预设消能减震体系的设计目标为：在罕遇地震下，安装阻尼器后结构最大层间位移角的减震率达20%以上，大部分楼层处于中等破坏水平，经过一般修理即可继续使用；在极罕遇地震下，结构不倒塌，结构最大层间位移角不超过1/50。根据《建筑地震破坏等级划分标准》，基本完好、轻微破坏、中等破坏、严重破坏、倒塌5个破坏等级对应的层间位移角限值分别为<1/550，(1.5-2)/550，(3-4)/550，<0.9/50，>1/50。

选取1999年Chi-Chi地震、1979年Imperial Valley地震和1971年San Fernando地震三条地震动时程对消能体系进行地震响应分析，三条地震动的反应谱平均值与规范谱在结构前三周期对应的最大误差为19.1%，均不超过规范要求的20%，满足规范中选波在统计意义上的要求，如图10所示。

图10 所选地震波与标准反应谱地震影响系数对比Fig.10 Comparison of seismic influence coefficient between selected seismic waves and standard response spectrum

4.2 罕遇地震作用下CRAD-VD消能体系地震响应

分别沿结构X，Y方向输入地震动，峰值加速度为0.4g。Imperial Valley波作用下层间位移角，如图11所示。三条地震波作用下无控结构、VD消能体系、CRAD-VD消能体系的层间位移角及减震率，如表2所示。

图11 Imperial Valley波(0.4g)地震作用下结构层间位移角Fig.11 Interlayer drift angle of the system under Imperial Valley wave (0.4g)

表2 罕遇地震作用下结构最大层间位移角及减震率Tab.2 Maximum inter-story drift angle and control effect of the structure under rarely occurred earthquake

由图11及表2可见，在X,Y方向，无控结构大部分楼层进入严重破坏状态，最大层间位移角达到1/66，接近倒塌；而安装VD和CRAD-VD后，结构层间位移角明显减小，大部分楼层处于中等破坏状态，减震率均达到20%以上；与VD消能体系相比，CRAD-VD消能体系达到并超过直接安装VD的控制效果，减震率均大于42%，说明由于CRAD的响应放大作用，串联小阻尼系数黏滞阻尼器的情况下，可获得优于直接安装大阻尼系数VD的位移控制效果。

Imperial Valley波作用下结构楼层位移和层间剪力,如图12所示。①安装VD和CRAD-VD后结构楼层位移均显著减少，说明阻尼器可有效耗散地震能量。不同地震波下，VD消能体系的减震率介于13%～39%；CRAD-VD消能体系减震率介于33%～60%，说明CRAD-VD对体系楼层位移的控制效果明显优于VD；②与无控结构相比，安装阻尼器VD和CRAD-VD后层间剪力均有明显的控制效果；与VD体系相比，CRAD-VD消能体系在局部楼层出现层间剪力略大的现象，但相对于无控的减震效果仍然较为显著。

在San Fernando波作用下，CRAD-VD与VD的累积滞回耗能时程曲线，如图13所示，总耗能和耗能比，如表3所示。可见，CRAD-VD的耗能能力是VD阻尼器的1.2～2.1倍，表明CRAD串联小阻尼系数的阻尼器具有比直接安装大阻尼系数阻尼器更大的耗能能力，体现了明显的力和耗能放大能力。

图12 Imperial Valley波(0.4g)地震作用下结构响应Fig.12 Seismic response under Imperial Valley wave (0.4g)

图13 San Fernando波作用下阻尼器累积滞回耗能时程曲线Fig.13 History curve of cumulative hysteretic energy dissipation of damper under San Fernando wave

表3 罕遇地震作用下阻尼器累积滞回耗能Tab.3 Cumulative hysteretic energy dissipation of damper under rarely occurred earthquake

在San Fernando波作用下，CRAD-VD、C-VD与VD的滞回曲线,如图14所示。C-VD与VD的位移时程曲线，如图15所示。由图14可见，X向和Y向阻尼器的位移幅值与阻尼力幅值均较接近，表明阻尼器参数设计及布置是合理的；与VD相比，CRAD串联阻尼器C-VD的阻尼力幅值小于VD阻尼力幅值，这是由于C-VD的阻尼系数仅为VD阻尼系数的0.4倍；但是，CRAD-VD的阻尼力幅值大于VD阻尼力幅值，滞回曲线呈现出较小位移幅值、较大阻尼力的现象，呈现出明显的响应放大效应。

图14 San Fernando波(0.4g)作用下阻尼器的滞回曲线Fig.14 Hysteretic curve of damper under San Fernando wave (0.4g)

图15 San Fernando波(0.4g)作用下VD及C-VD的位移时程Fig.15 Displacement time history of VD and C-VD under San Fernando wave (0.4g)

由图14和图15可见，CRAD串联阻尼器C-VD始终在设定凸轮偏心距e=0.02 m内做往复运动。

无控装配式结构在罕遇地震作用下，大部分的框架柱及框架梁都出现了严重损伤，节点破坏数量达到了867个，结构面临倒塌的风险；VD体系和CRAD-VD体系框架柱大部分处于中等破坏状态，节点破坏数量分别为214个和189个，说明安装响应阻尼器后结构的安全性大幅度提高，有利于震后救灾的实现；与直接安装大吨位VD相比，CRAD-VD仅需串联小吨位C-VD就可达到大吨位VD类似的效果，能有效降低框架梁柱在罕遇地震下的损伤程度。

4.3 极罕遇地震下CRAD-VD消能体系响应

在Imperial Valley极罕遇地震动(0.64g)作用下，无控结构、VD消能体系、CRAD-VD消能体系的层间位移角，如图16所示，三条地震波作用下的层间位移角及减震率，如表4所示。

图16 Imperial Valley波(0.64g)作用下结构层间位移角Fig.16 Interlayer drift angle of the system under Imperial Valley wave (0.64g)

表4 极罕遇地震作用下结构最大层间位移角及减震率Tab.4 Maximum inter-story drift angle and control effect of the structure under extremely rare occurred earthquake

由图16和表4可见，无控结构均处于严重破坏甚至倒塌状态。Imperial Valley(Y向)和San Fernando波作用下，VD体系处于严重破坏状态。而CRAD-VD体系处于中等破坏状态，显著超过VD体系的控制效果，说明CRAD在极罕遇地震下仍然能有效发挥其阻尼放大效应。此外，极罕遇地震作用下阻尼器累积滞回耗能和耗能比，如表5所示。可见，CRAD-VD的耗能能力是VD的1.35～1.87倍，有效地降低了框架梁柱的损伤程度，说明串联小阻尼系数阻尼器的CRAD-VD耗散了更多的地震能量，提高了装配整体式混凝土框架结构在极罕遇地震作用下的抗震安全性。

San Fernando波作用下阻尼器的滞回曲线，如图17所示。可见，极罕遇地震作用下，VD消能体系中阻尼器的位移远超20 mm，而CRAD-VD所串联阻尼器的位移在往复运动过程中始终小于或等于预设偏心距20 mm，有效地避免了阻尼器位移超限而失效的问题。

图17 San Fernando波(0.64g)作用下串联阻尼器的滞回曲线Fig.17 Hysteretic curve of C-VD under San Fernando wave (0.64g)

表5 极罕遇地震作用下阻尼器累积滞回耗能Tab.5 Cumulative hysteretic energy dissipation of damper under extremely rare occurred earthquake kJ

5 结 论

为便于采用ABAQUS软件对CRAD-VD消能减震三维结构进行抗震性能分析，基于ABAQUS开发了CRAD-VD耗能单元，编制了CRAD-VD消能体系非线性分析的显式求解算法，并与MATLAB计算程序进行了对比验证，进而开展了CRAD-VD消能三维减震体系在罕遇和极罕遇地震作用下的抗震性能分析，结果表明：

(1)CRAD-VD滞回曲线呈现高度非线性，采用隐式算法进行动力分析可能会出现由于迭代次数过多、未达到收敛要求而终止计算的情况，显式算法无需迭代，只要增量步长小于稳定极限值即不存在收敛问题，更适合CRAD-VD消能减震结构的计算分析。基于ABAQUS编制的显式算法计算消能体系的反应结果均与相应的MATLAB程序结果吻合，证明了二次开发的准确性和可靠性。

(2)在罕遇和极罕遇地震作用下，无控装配式框架结构、单独设置阻尼器VD减震体系、安装CRAD-VD减震体系的响应结果表明，CRAD-VD串联小阻尼系数的黏滞阻尼器取得了相似甚至优于安装大阻尼系数阻尼器的减震效果，且有效地降低了装配式框架梁柱的损伤程度，确保了结构的抗震安全，充分体现了CRAD的响应放大效应。

(3)极罕遇地震作用下，CRAD-VD串联阻尼器的位移幅值在往复运动过程中始终等于预设偏心距，有效地避免了阻尼器位移超限而失效的问题。