魯 正，呂西林

（1.同濟大學 土木工程防災國家重點實驗室，上海200092；2.同濟大學 結(jié)構(gòu)工程與防災研究所，上海200092）

顆粒阻尼器將裝有金屬或其他材料等顆粒的容器附著在結(jié)構(gòu)振動較大的部位，利用顆粒之間以及顆粒與容器壁之間的非彈性碰撞和摩擦來消耗系統(tǒng)振動能量，是一種附加質(zhì)量式被動阻尼器.該技術(shù)具有概念簡單、減振頻帶寬、溫度不敏感、耐久性好、易于用在惡劣環(huán)境等優(yōu)點，具有在土木工程應用的良好前景.Naeim 等［1］介紹了一幢位于圣地亞哥市中心經(jīng)受了2010年智利地震考驗的附加顆粒阻尼器的高層建筑，該建筑在地震后結(jié)構(gòu)完好，沒有任何破壞；魯正，呂西林等［2］完成了附加顆粒阻尼器的三層鋼框架的振動臺試驗，表明附加很小質(zhì)量比的顆粒阻尼器能夠有效減小主體結(jié)構(gòu)的響應，尤其是均方根響應.

考慮到顆粒碰撞帶來的非線性行為，系統(tǒng)很難求得解析解，因而學者們研究出一些簡化方法和數(shù)值方法.Papalou和Masri［3］把多顆粒阻尼器簡化等效為等質(zhì)量的單顆粒阻尼器；Friend和Kinra［4］通過把多顆粒模擬為一個凝聚的質(zhì)量塊，把各種機理引起的能量耗散打包在一個新的概念“有效恢復系數(shù)”，該系數(shù)通過實驗擬合得到，從而提出一種解析方法.此外，運用于顆粒阻尼模擬的方法還有：回歸模型法［5］，恢復力曲面法［6］和功率輸入法［7］等.盡管這些簡化模型和基于實驗的研究取得了不小的成果，然而這些終究是基于現(xiàn)象的方法，結(jié)論也很難推廣到除了該實驗之外的其他情形.近年來，離散單元法（discrete element method，DEM）被引入到顆粒阻尼器的分析.由于該方法能模擬顆粒之間，以及顆粒與容器壁之間的相互作用，因此能更合理地定量分析顆粒阻尼器的性能.本文介紹基于離散單元法的顆粒阻尼器對多自由度結(jié)構(gòu)進行減震控制的數(shù)值模擬方法，并進行試驗驗證.

1 數(shù)值模擬方法

離散單元法［8］是一種按時步迭代求解，研究離散體力學行為的數(shù)值分析方法，它把離散體劃分為眾多離散單元的集合，根據(jù)接觸定律和牛頓第二定律描述其運動.該方法認為只要時步取值足夠小，則在該時步內(nèi)，單元的擾動只會傳播到與其相鄰的單元，而不足以傳播到其他更遠的單元.據(jù)此，作用在某一單元上的外力就可以通過與其相鄰單元的相互作用情況求得，進而求得整個離散體的整體運動形態(tài).

1.1 控制方程

圖1a示意了一個在頂層附加顆粒阻尼器的多自由度結(jié)構(gòu)，N為樓層數(shù)，其控制方程為

式中：M，C，K分別為質(zhì)量，阻尼和剛度矩陣分別為位移，速度和加速度矩陣；E為慣性質(zhì)量矩陣為地面加速度；F為顆粒對結(jié)構(gòu)的接觸力向量，這也是顆粒與主體結(jié)構(gòu)之間的聯(lián)系紐帶.

對于顆粒i，某一時刻的控制方程為

式中：mi為顆粒的質(zhì)量；Ii為顆粒的慣性矩；g為 重力加速度向量；pi，φi分別為顆粒的位置向量和角位移向量分別為顆粒的平動加速度向量和角加速度向量為顆粒i和顆粒j之間的法向接觸力和切向接觸力（若顆粒i與容器壁接觸，則j代表容器壁）.接觸力作用在兩個顆粒的接觸點，而不是在顆粒的質(zhì)心，切向接觸力會產(chǎn)生扭矩Tij，使顆粒產(chǎn)生旋轉(zhuǎn).對半徑為r的球形顆粒其中，nij為顆粒i的質(zhì)心指向顆粒j的質(zhì)心的單位向量，“×”表示向量叉積，ki為與顆粒i相接觸的顆粒數(shù)目.

1.2 接觸力模型

接觸力模型用來定量確定法向力和切向力.本文采用計算速度快且概念簡單的線性模型，即法向為線性接觸力模型，切向為庫倫摩擦力模型，作為對接觸力模型選取的有益探討.

圖1b是顆粒與容器壁的法向線性接觸力模型.圖中，k2為彈簧 剛度為角頻率，通過合理選擇ω2的值來模擬容器壁，Masri采用圖1b的接觸力模型，對單自由度體系附加沖擊阻尼器的系統(tǒng)進行了解析解和數(shù)值解的對比分析［9］，指出通過調(diào)整彈簧剛度ω2／ωn≥20 來模擬容器壁是合適的.c2為阻尼系數(shù)，ζ2＝c2／2mω2，為阻尼比，能用來模擬非彈性碰撞，因此各種恢復系數(shù)（coefficient of restitution，兩物體碰撞后的相對速度和碰撞前的相對速度的比值的絕對值）可以通過調(diào)整ζ2 的值來實現(xiàn).顆粒之間的法向線性接觸力模型也類似，用ω3，c3和ζ3 代表顆粒間模擬法向彈簧剛度的角頻率，法向阻尼單元的阻尼系數(shù)和阻尼比.從而，法向力表示為

式中：δn和是顆粒i相對于顆粒j的位移和速度；是顆粒與容器壁的距離.采用庫侖摩擦力模型，切向接觸力表示為

圖1 計算模型簡圖Fig.1 Schematic diagram of computational models

式中：μs 是顆粒間或者顆粒與容器壁之間的摩擦系數(shù)是顆粒i相對于顆粒j的切向速度.

1.3 法向阻尼系數(shù)

法向阻尼系數(shù)可以由球和容器壁碰撞的恢復系數(shù)導出，碰撞模型見圖1b.球與容器壁接觸后的動力方程為

顆粒的初始位置為x0＝0，初始速度為，接觸前的入射速度為，求解得到位移響應和速度響應為

根據(jù)顆粒恢復系數(shù)e的定義：從而可以得到法向阻尼系數(shù)和顆?；謴拖禂?shù)的關(guān)系，通過調(diào)整ζ2 的值就可以模擬各種材料顆粒的彈性狀態(tài)，如圖2所示.

圖2 法向阻尼系數(shù)與恢復系數(shù)關(guān)系圖Fig.2 Relationship between damping ratio and coefficient of restitution

1.4 計算時步

在離散元的數(shù)值模擬中，計算時步對整個計算結(jié)果的穩(wěn)定性具有十分重要的作用，若時步取值太大，則難以捕捉到顆粒碰撞的力學行為；若時步取值太小，則計算時間大大增長.當計算時步大約為碰撞接觸時間的1／5時，該碰撞的力學特性能被合理模擬并保證計算的穩(wěn)定性.在典型的離散元計算中，實際計算時步通常取為臨界時步的一定比例.臨界時步表示為［10］

式中：kmax是單元間的最大剛度；mmin是單元的最小質(zhì)量，于是是系統(tǒng)的最大自振頻率；α是用戶自定義參數(shù)，計算經(jīng)驗表明該值在0.1時所得到的計算時步可以保證計算穩(wěn)定［11］.

1.5 接觸檢測算法

離散元計算的過程中，在每一個計算時步，都要進行接觸判斷，因此接觸檢測算法的合理與否直接關(guān)系到離散元模擬的效率.對于n個顆粒，最直觀的接觸判斷方法是遍歷判斷，即每一個顆粒均與除自己之外的其他顆粒做對比，以此來判斷是否接觸，這種判斷算法的復雜度為n2.一般離散元法常用的檢測算法是采用“盒式”判斷，即將整個空間分割成多個立方體式的小空間，在每個小空間內(nèi)，一個或者更多個顆粒相互運動，在進行接觸判斷時，只需要判斷某個顆粒與該小空間內(nèi)的其他顆粒的相對位置即可，而不需要判別此顆粒與其他小空間內(nèi)的顆粒的相對位置，該種判斷的復雜度為nln（n）.本研究采用一種更為高效的由Munjiza提出的非二元搜索算法（no binary search，NBS）［12］，該方法將整個空間劃分成邊長為最小顆粒直徑的若干個正方體，以保證每個顆粒可以且僅可以位于一個正方體之內(nèi)，該方法的復雜度僅為n，最適用于顆粒半徑相差不大的情況，而本研究所采用的顆粒阻尼器為大小一致的顆粒，因而采用NBS方法是最合理有效的.

1.6 數(shù)值模擬流程

應用離散單元法模擬附加顆粒阻尼器的多自由度結(jié)構(gòu)的過程簡述如下：首先，判斷顆粒之間，顆粒與容器壁之間的相對位置，若δn＞0，作用在顆粒上的接觸力可以通過式（3）和（4）求得，若δn≤0，則無接觸力；其次，對作用在一個顆粒上的所有的接觸力求和，包括顆粒之間的接觸力和顆粒與容器壁的接觸力；再次，顆粒的運動可以通過式（2）求得；以上過程對所有顆粒順次進行；最后，對所有顆粒與容器壁的接觸力累加，就得到等式（1）中的力F，對其求解，即得到主體結(jié)構(gòu)的響應.

2 試驗驗證

為了驗證本文提出的基于離散元的顆粒阻尼器減震控制模擬方法的可行性和合理性，開展了一個三層鋼框架附加顆粒阻尼器的模擬地震振動臺試驗.

2.1 試驗簡介

試驗采用的結(jié)構(gòu)模型是一個三層鋼框架，并在頂層附加顆粒阻尼器.為了模擬土木工程中的高層建筑，調(diào)節(jié)框架的基頻至1.0 Hz左右，在結(jié)構(gòu)各層附加質(zhì)量塊.試驗時，第一至第三層的實際質(zhì)量（包含結(jié)構(gòu)自重）分別為1 915，1 915和2 124kg.該框架的阻尼比ζ1 為0.013，前三階自振頻率分別為1.07，3.2和4.8Hz.顆粒阻尼器由鋼板組成，分為4個相同并且沿著震動方向?qū)ΨQ的立方體鐵盒，長×寬×高尺寸為：0.49m×0.49m×0.5m.每個鐵盒子內(nèi)分別放置63個鋼球，直徑為50.8 mm，總計質(zhì)量135kg，占系統(tǒng)總體質(zhì)量的2.25%.關(guān)于試驗過程和結(jié)果在文獻［2］中有詳細介紹，本文僅介紹其數(shù)值模擬參數(shù)和結(jié)果.

根據(jù)試驗記錄，從采集的臺面加速度時程曲線中截取一個完整波作為計算輸入波，將依據(jù)此波計算出來的結(jié)果與同時采集的模型頂層位移時程曲線相比較，通過理論與實測曲線的符合程度來驗證所提的數(shù)值模型和算法的正確性.在進行數(shù)值模擬時，為方便起見，將各質(zhì)點的位移和速度初值賦為零.因此，計算結(jié)果與試驗記錄相比較，在起始段會有一定出入.隨著計算過程的逐步進行，初值選取對結(jié)構(gòu)體系反應的影響將逐步減小.

2.2 參數(shù)取值

2.2.1 摩擦系數(shù)

本文主要參考Wong等［13］在研究顆粒阻尼器的運動時，特地開展的測定鋼球摩擦系數(shù)的試驗，該試驗對不同速度不同直徑的鋼球與鋼板之間的摩擦系數(shù)進行測定，在球速較快的狀況下，摩擦系數(shù)均值為0.47.由于顆粒之間接觸區(qū)域很小，可以看作為點接觸，同樣地，顆粒與容器壁之間的接觸也可以認為是點接觸，因此，當顆粒與容器的材料相同，表面狀況相同時，顆粒間的摩擦系數(shù)可以看作與顆?！萜鞅陂g的摩擦系數(shù)相同或者非常相近，本文近似取為相等［13］.

2.2.2 模擬碰撞的法向阻尼單元的阻尼比

根據(jù)葛藤等［14］的研究，在鋼球速度為0.1～1m·s－1時，對應的恢復系數(shù)為0.7～0.8，這與Xia等［15］在模擬中取用恢復系數(shù)為0.78吻合.由恢復系數(shù)與法向阻尼單元的阻尼系數(shù)的對應關(guān)系圖2 得，對應的阻尼比為0.1.

2.2.3 模擬碰撞的法向彈簧剛度和計算時步

如1.2節(jié)討論，法向彈簧剛度的取值為105N·m－1，滿足ω2／ωn≥20的要求.相應地，當計算時步取為1×10－4時，也滿足1.4 節(jié)對計算時步取值的要求，計算得以穩(wěn)定.

2.2.4 試驗框架特性

未附加阻尼器的空框架的質(zhì)量矩陣通過每一層的實際質(zhì)量來建立.剛度矩陣首先根據(jù)結(jié)構(gòu)力學公式，計算每一層所有柱子的抗側(cè)剛度來得到.此外，在試驗前，通過白噪聲掃頻得到結(jié)構(gòu)的自振特性和阻尼比，考慮到模型加工誤差，抗側(cè)剛度計算公式的理想假定與實際制作的誤差等因素，最終根據(jù)前三階框架自振頻率對結(jié)構(gòu)的剛度矩陣進行適當修正，得到最終用于計算的剛度矩陣.阻尼矩陣根據(jù)掃頻得到的阻尼比求得.

綜上，試驗框架的質(zhì)量M，kg；、剛度K，N·m－1；阻尼比ζ1 如下，用于計算的系統(tǒng)參數(shù)見表1.

表1 系統(tǒng)參數(shù)取值Tab.1 Values of system parameters

2.3 數(shù)值模擬結(jié)果

圖3畫出了附加顆粒阻尼器的框架模型頂層在0.2g（g為重力加速度）地震激勵（Kobe波）下的位移和加速度響應的計算值和試驗值的對比曲線，可以看到兩者符合較好.其中，位移時程的符合度更好，這是因為位移是加速度的積分，因此其曲線更光滑.對于加速度曲線，可以看到在顆粒與主體結(jié)構(gòu)碰撞的時候，會有突變.表2列出了模型頂層峰值位移在不同地震激勵下的計算值與試驗值的對比，發(fā)現(xiàn)兩者也吻合良好；此外，還列出了位移時程曲線誤差序列的均方根值，在圖3曲線直觀對比的基礎上，進一步提供一些定量的概念；考慮到附加顆粒阻尼器的多自由度結(jié)構(gòu)在地震作用下，是一個復雜的非線性系統(tǒng)；而且數(shù)值模擬時，為方便起見，將各質(zhì)點的位移和速度初值賦為零，導致在計算結(jié)果與試驗記錄相比的初始段會有一定出入，但是隨著計算過程的逐步進行，初值選取對結(jié)構(gòu)體系反應的影響（尤其是峰值響應的影響）逐步減小等因素，可以認為這些偏差是可以接受的.這些都說明本文提出的數(shù)值模擬方法能夠較合理地計算出附加顆粒阻尼器系統(tǒng)在實際地震激勵下的響應.當然了，對于該系統(tǒng)的更精確有效的模擬，還需要進一步的研究，尤其是對接觸模型以及參數(shù)取值的研究.

圖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 結(jié) 論

本文基于離散單元法，建立了顆粒阻尼器對多自由度結(jié)構(gòu)進行減震控制的數(shù)值模擬方法，并介紹了系統(tǒng)參數(shù)的確定方法.該方法采用合理的計算時步和高效的接觸檢測算法，能較好地模擬顆粒間以及顆粒與容器壁之間的相互作用，包括摩擦力，重力，碰撞力等.為了驗證該方法的可行性和準確性，開展了附加顆粒阻尼器的多層框架的振動臺試驗，通過試驗結(jié)果與數(shù)值模擬結(jié)果的對比發(fā)現(xiàn)兩者吻合良好，說明該數(shù)值方法可以較好地模擬顆粒阻尼器的減震控制效果，為進一步了解顆粒阻尼器的物理本質(zhì)，并進而推廣其在土木工程中的應用，提供了一條可行且高效的數(shù)值分析途徑.

［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］ 劉雁梅，黃協(xié)清，陳天寧.非阻塞性顆粒阻尼加筋板振動功率流的研究［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］ 葛藤，賈智宏，周克棟.鋼球和剛性平面彈塑性正碰撞恢復系數(shù)研究［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.