颗粒阻尼器减震控制的数值模拟

鲁 正，吕西林

（1.同济大学 土木工程防灾国家重点实验室，上海200092；2.同济大学 结构工程与防灾研究所，上海200092）

颗粒阻尼器将装有金属或其他材料等颗粒的容器附着在结构振动较大的部位，利用颗粒之间以及颗粒与容器壁之间的非弹性碰撞和摩擦来消耗系统振动能量，是一种附加质量式被动阻尼器.该技术具有概念简单、减振频带宽、温度不敏感、耐久性好、易于用在恶劣环境等优点，具有在土木工程应用的良好前景.Naeim 等［1］介绍了一幢位于圣地亚哥市中心经受了2010年智利地震考验的附加颗粒阻尼器的高层建筑，该建筑在地震后结构完好，没有任何破坏；鲁正，吕西林等［2］完成了附加颗粒阻尼器的三层钢框架的振动台试验，表明附加很小质量比的颗粒阻尼器能够有效减小主体结构的响应，尤其是均方根响应.

考虑到颗粒碰撞带来的非线性行为，系统很难求得解析解，因而学者们研究出一些简化方法和数值方法.Papalou和Masri［3］把多颗粒阻尼器简化等效为等质量的单颗粒阻尼器；Friend和Kinra［4］通过把多颗粒模拟为一个凝聚的质量块，把各种机理引起的能量耗散打包在一个新的概念“有效恢复系数”，该系数通过实验拟合得到，从而提出一种解析方法.此外，运用于颗粒阻尼模拟的方法还有：回归模型法［5］，恢复力曲面法［6］和功率输入法［7］等.尽管这些简化模型和基于实验的研究取得了不小的成果，然而这些终究是基于现象的方法，结论也很难推广到除了该实验之外的其他情形.近年来，离散单元法（discrete element method，DEM）被引入到颗粒阻尼器的分析.由于该方法能模拟颗粒之间，以及颗粒与容器壁之间的相互作用，因此能更合理地定量分析颗粒阻尼器的性能.本文介绍基于离散单元法的颗粒阻尼器对多自由度结构进行减震控制的数值模拟方法，并进行试验验证.

1 数值模拟方法

离散单元法［8］是一种按时步迭代求解，研究离散体力学行为的数值分析方法，它把离散体划分为众多离散单元的集合，根据接触定律和牛顿第二定律描述其运动.该方法认为只要时步取值足够小，则在该时步内，单元的扰动只会传播到与其相邻的单元，而不足以传播到其他更远的单元.据此，作用在某一单元上的外力就可以通过与其相邻单元的相互作用情况求得，进而求得整个离散体的整体运动形态.

1.1 控制方程

图1a示意了一个在顶层附加颗粒阻尼器的多自由度结构，N为楼层数，其控制方程为

式中：M，C，K分别为质量，阻尼和刚度矩阵分别为位移，速度和加速度矩阵；E为惯性质量矩阵为地面加速度；F为颗粒对结构的接触力向量，这也是颗粒与主体结构之间的联系纽带.

对于颗粒i，某一时刻的控制方程为

式中：mi为颗粒的质量；Ii为颗粒的惯性矩；g为 重力加速度向量；pi，φi分别为颗粒的位置向量和角位移向量分别为颗粒的平动加速度向量和角加速度向量为颗粒i和颗粒j之间的法向接触力和切向接触力（若颗粒i与容器壁接触，则j代表容器壁）.接触力作用在两个颗粒的接触点，而不是在颗粒的质心，切向接触力会产生扭矩Tij，使颗粒产生旋转.对半径为r的球形颗粒其中，nij为颗粒i的质心指向颗粒j的质心的单位向量，“×”表示向量叉积，ki为与颗粒i相接触的颗粒数目.

1.2 接触力模型

接触力模型用来定量确定法向力和切向力.本文采用计算速度快且概念简单的线性模型，即法向为线性接触力模型，切向为库伦摩擦力模型，作为对接触力模型选取的有益探讨.

图1b是颗粒与容器壁的法向线性接触力模型.图中，k2为弹簧 刚度为角频率，通过合理选择ω2的值来模拟容器壁，Masri采用图1b的接触力模型，对单自由度体系附加冲击阻尼器的系统进行了解析解和数值解的对比分析［9］，指出通过调整弹簧刚度ω2／ωn≥20 来模拟容器壁是合适的.c2为阻尼系数，ζ2＝c2／2mω2，为阻尼比，能用来模拟非弹性碰撞，因此各种恢复系数（coefficient of restitution，两物体碰撞后的相对速度和碰撞前的相对速度的比值的绝对值）可以通过调整ζ2 的值来实现.颗粒之间的法向线性接触力模型也类似，用ω3，c3和ζ3 代表颗粒间模拟法向弹簧刚度的角频率，法向阻尼单元的阻尼系数和阻尼比.从而，法向力表示为

式中：δn和是颗粒i相对于颗粒j的位移和速度；是颗粒与容器壁的距离.采用库仑摩擦力模型，切向接触力表示为

图1 计算模型简图Fig.1 Schematic diagram of computational models

式中：μs 是颗粒间或者颗粒与容器壁之间的摩擦系数是颗粒i相对于颗粒j的切向速度.

1.3 法向阻尼系数

法向阻尼系数可以由球和容器壁碰撞的恢复系数导出，碰撞模型见图1b.球与容器壁接触后的动力方程为

颗粒的初始位置为x0＝0，初始速度为，接触前的入射速度为，求解得到位移响应和速度响应为

根据颗粒恢复系数e的定义：从而可以得到法向阻尼系数和颗粒恢复系数的关系，通过调整ζ2 的值就可以模拟各种材料颗粒的弹性状态，如图2所示.

图2 法向阻尼系数与恢复系数关系图Fig.2 Relationship between damping ratio and coefficient of restitution

1.4 计算时步

在离散元的数值模拟中，计算时步对整个计算结果的稳定性具有十分重要的作用，若时步取值太大，则难以捕捉到颗粒碰撞的力学行为；若时步取值太小，则计算时间大大增长.当计算时步大约为碰撞接触时间的1／5时，该碰撞的力学特性能被合理模拟并保证计算的稳定性.在典型的离散元计算中，实际计算时步通常取为临界时步的一定比例.临界时步表示为［10］

式中：kmax是单元间的最大刚度；mmin是单元的最小质量，于是是系统的最大自振频率；α是用户自定义参数，计算经验表明该值在0.1时所得到的计算时步可以保证计算稳定［11］.

1.5 接触检测算法

离散元计算的过程中，在每一个计算时步，都要进行接触判断，因此接触检测算法的合理与否直接关系到离散元模拟的效率.对于n个颗粒，最直观的接触判断方法是遍历判断，即每一个颗粒均与除自己之外的其他颗粒做对比，以此来判断是否接触，这种判断算法的复杂度为n2.一般离散元法常用的检测算法是采用“盒式”判断，即将整个空间分割成多个立方体式的小空间，在每个小空间内，一个或者更多个颗粒相互运动，在进行接触判断时，只需要判断某个颗粒与该小空间内的其他颗粒的相对位置即可，而不需要判别此颗粒与其他小空间内的颗粒的相对位置，该种判断的复杂度为nln（n）.本研究采用一种更为高效的由Munjiza提出的非二元搜索算法（no binary search，NBS）［12］，该方法将整个空间划分成边长为最小颗粒直径的若干个正方体，以保证每个颗粒可以且仅可以位于一个正方体之内，该方法的复杂度仅为n，最适用于颗粒半径相差不大的情况，而本研究所采用的颗粒阻尼器为大小一致的颗粒，因而采用NBS方法是最合理有效的.

1.6 数值模拟流程

应用离散单元法模拟附加颗粒阻尼器的多自由度结构的过程简述如下：首先，判断颗粒之间，颗粒与容器壁之间的相对位置，若δn＞0，作用在颗粒上的接触力可以通过式（3）和（4）求得，若δn≤0，则无接触力；其次，对作用在一个颗粒上的所有的接触力求和，包括颗粒之间的接触力和颗粒与容器壁的接触力；再次，颗粒的运动可以通过式（2）求得；以上过程对所有颗粒顺次进行；最后，对所有颗粒与容器壁的接触力累加，就得到等式（1）中的力F，对其求解，即得到主体结构的响应.

2 试验验证

为了验证本文提出的基于离散元的颗粒阻尼器减震控制模拟方法的可行性和合理性，开展了一个三层钢框架附加颗粒阻尼器的模拟地震振动台试验.

2.1 试验简介

试验采用的结构模型是一个三层钢框架，并在顶层附加颗粒阻尼器.为了模拟土木工程中的高层建筑，调节框架的基频至1.0 Hz左右，在结构各层附加质量块.试验时，第一至第三层的实际质量（包含结构自重）分别为1 915，1 915和2 124kg.该框架的阻尼比ζ1 为0.013，前三阶自振频率分别为1.07，3.2和4.8Hz.颗粒阻尼器由钢板组成，分为4个相同并且沿着震动方向对称的立方体铁盒，长×宽×高尺寸为：0.49m×0.49m×0.5m.每个铁盒子内分别放置63个钢球，直径为50.8 mm，总计质量135kg，占系统总体质量的2.25%.关于试验过程和结果在文献［2］中有详细介绍，本文仅介绍其数值模拟参数和结果.

根据试验记录，从采集的台面加速度时程曲线中截取一个完整波作为计算输入波，将依据此波计算出来的结果与同时采集的模型顶层位移时程曲线相比较，通过理论与实测曲线的符合程度来验证所提的数值模型和算法的正确性.在进行数值模拟时，为方便起见，将各质点的位移和速度初值赋为零.因此，计算结果与试验记录相比较，在起始段会有一定出入.随着计算过程的逐步进行，初值选取对结构体系反应的影响将逐步减小.

2.2 参数取值

2.2.1 摩擦系数

本文主要参考Wong等［13］在研究颗粒阻尼器的运动时，特地开展的测定钢球摩擦系数的试验，该试验对不同速度不同直径的钢球与钢板之间的摩擦系数进行测定，在球速较快的状况下，摩擦系数均值为0.47.由于颗粒之间接触区域很小，可以看作为点接触，同样地，颗粒与容器壁之间的接触也可以认为是点接触，因此，当颗粒与容器的材料相同，表面状况相同时，颗粒间的摩擦系数可以看作与颗粒—容器壁间的摩擦系数相同或者非常相近，本文近似取为相等［13］.

2.2.2 模拟碰撞的法向阻尼单元的阻尼比

根据葛藤等［14］的研究，在钢球速度为0.1～1m·s－1时，对应的恢复系数为0.7～0.8，这与Xia等［15］在模拟中取用恢复系数为0.78吻合.由恢复系数与法向阻尼单元的阻尼系数的对应关系图2 得，对应的阻尼比为0.1.

2.2.3 模拟碰撞的法向弹簧刚度和计算时步

如1.2节讨论，法向弹簧刚度的取值为105N·m－1，满足ω2／ωn≥20的要求.相应地，当计算时步取为1×10－4时，也满足1.4 节对计算时步取值的要求，计算得以稳定.

2.2.4 试验框架特性

未附加阻尼器的空框架的质量矩阵通过每一层的实际质量来建立.刚度矩阵首先根据结构力学公式，计算每一层所有柱子的抗侧刚度来得到.此外，在试验前，通过白噪声扫频得到结构的自振特性和阻尼比，考虑到模型加工误差，抗侧刚度计算公式的理想假定与实际制作的误差等因素，最终根据前三阶框架自振频率对结构的刚度矩阵进行适当修正，得到最终用于计算的刚度矩阵.阻尼矩阵根据扫频得到的阻尼比求得.

综上，试验框架的质量M，kg；、刚度K，N·m－1；阻尼比ζ1 如下，用于计算的系统参数见表1.

表1 系统参数取值Tab.1 Values of system parameters

2.3 数值模拟结果

图3画出了附加颗粒阻尼器的框架模型顶层在0.2g（g为重力加速度）地震激励（Kobe波）下的位移和加速度响应的计算值和试验值的对比曲线，可以看到两者符合较好.其中，位移时程的符合度更好，这是因为位移是加速度的积分，因此其曲线更光滑.对于加速度曲线，可以看到在颗粒与主体结构碰撞的时候，会有突变.表2列出了模型顶层峰值位移在不同地震激励下的计算值与试验值的对比，发现两者也吻合良好；此外，还列出了位移时程曲线误差序列的均方根值，在图3曲线直观对比的基础上，进一步提供一些定量的概念；考虑到附加颗粒阻尼器的多自由度结构在地震作用下，是一个复杂的非线性系统；而且数值模拟时，为方便起见，将各质点的位移和速度初值赋为零，导致在计算结果与试验记录相比的初始段会有一定出入，但是随着计算过程的逐步进行，初值选取对结构体系反应的影响（尤其是峰值响应的影响）逐步减小等因素，可以认为这些偏差是可以接受的.这些都说明本文提出的数值模拟方法能够较合理地计算出附加颗粒阻尼器系统在实际地震激励下的响应.当然了，对于该系统的更精确有效的模拟，还需要进一步的研究，尤其是对接触模型以及参数取值的研究.

图3 附加颗粒阻尼器的框架模型顶层在0.2g 地震激励下的响应（Kobe波）Fig.3 Response at roof level of the test frame with particle dampers under 0.2g earthquake excitations（Kobe Wave）

表2 模型顶层位移响应计算值与试验值对比Tab.2 Comparison of the calculated and experimental results for the displacement on the roof of the test frame

3 结 论

本文基于离散单元法，建立了颗粒阻尼器对多自由度结构进行减震控制的数值模拟方法，并介绍了系统参数的确定方法.该方法采用合理的计算时步和高效的接触检测算法，能较好地模拟颗粒间以及颗粒与容器壁之间的相互作用，包括摩擦力，重力，碰撞力等.为了验证该方法的可行性和准确性，开展了附加颗粒阻尼器的多层框架的振动台试验，通过试验结果与数值模拟结果的对比发现两者吻合良好，说明该数值方法可以较好地模拟颗粒阻尼器的减震控制效果，为进一步了解颗粒阻尼器的物理本质，并进而推广其在土木工程中的应用，提供了一条可行且高效的数值分析途径.

［1］ Naeim F，Lew M，Carpenter L D，et al.Performance of tall buildings in Santiago，Chile during the 27 February 2010 offshore Maule，Chile earthquake［J］.The Structural Design of Tall and Special Buildings，2011，20（1）：1.

［2］ 鲁正，吕西林，闫维明.颗粒阻尼器减震控制的试验研究［J］.土木工程学报，2012，45（增刊1）：243.LU Zheng， LU Xilin， YAN Weiming. Experimental investigation into the vibration control effects of particle dampers［J］.China Civil Engineering Journal，2012，45（Supplement 1）：243.

［3］ Papalou A，Masri S F.Response of impact dampers with granular materials under random excitation ［J］.Earthquake Engineering ＆Structural Dynamics，1996，25（3）：253.

［4］ Friend R D，Kinra V K.Particle impacting damping ［J］.Journal of Sound and Vibration，2000，233（1）：93.

［5］ 胡溧，黄其柏，许智生.颗粒阻尼的回归分析研究［J］.中国机械工程，2008，19（23）：2834.HU Li，HUANG Qibai，XU Zhisheng.Regression analysis of particle damping［J］.Chinese Mechanical Engineering，2008，19（23）：2834.

［6］ 蒋华，陈前.恢复力曲面法在颗粒阻尼器研究中的应用［J］.振动、测试与诊断，2007，27（3）：228.JIANG Hua，CHEN Qian. Application of resoting force superface method to particle damping research［J］.Journal of Vibration，Measurement ＆Diagnosis，2007，27（3）：228.

［7］ 刘雁梅，黄协清，陈天宁.非阻塞性颗粒阻尼加筋板振动功率流的研究［J］.西安交通大学学报，2001，35（1）：61.LIU Yanmei，HUANG Xieqing，CHEN Tianning.Study on power flow of a stiffened plate with non－obstructive particle damping［J］.Journal of Xi’an Jiaotong University，2001，35（1）：61.

［8］ Cundall P A，Strack O.A distinct element model for granular assemblies［J］.Geotechnique，1979，29：47.

［9］ Masri S F，Ibrahim A M.Stochastic excitation of a simple system with impact damper［J］.Earthquake Engineering ＆Structural Dynamics，1973，1（4）：337.

［10］ Tavarez F A.Discrete element method for modeling solid and particulate materials［D］.Madison：Engineering Mechanics of University of Wisconsin－Madison，2005.

［11］ Jensen R P，Bosscher P J，Plesha M E.DEM simulation of granular media－structure interface： effects of surface roughness and particle shape ［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1999（23）：531.

［12］ Munjiza A，Andrews K R F.NBS contact detection algorithm for bodies of similar size ［J］.International Journal for Numerical Methods in Engineering，1998，43：131.

［13］ Wong C X，Daniel M C，Rongong J A.Energy dissipation prediction of particle dampers ［J］.Journal of Sound and Vibration，2009，319（1／2）：91.

［14］ 葛藤，贾智宏，周克栋.钢球和刚性平面弹塑性正碰撞恢复系数研究［J］.工程力学，2008，25（6）：209.GE Teng，JIA Zhihong， ZHOU Kedong. Research on elastoplastic normal impact of steel spheres against a rigid plane［J］.Engineering Mechanics，2008，25（6）：209.

［15］ Xia Z W，Liu X D，Shan Y C，et al.Coupling simulation algorithm of discrete element method and finite element method for particle damper［J］.Journal of Low Frequency Noise，Vibration and Active Control，2009，28（3）：197.