顆粒物質流變學行為和材料參數對顆粒阻尼器能量耗散的影響

蘇 凡， 張 航， 尹忠俊

(北京科技大學 機械工程學院， 北京 100083)

顆粒阻尼技術，作為一種有效的振動控制手段，被廣泛地應用于航空航天、機械以及土木等領域[1]。顆粒阻尼的原理是在振動結構體上加工一定數量的孔洞，在其中填充適當數量的金屬、非金屬顆粒或通過將裝填有顆粒的空腔附著在結構振動較大的區域，利用顆粒與顆粒以及顆粒與空腔之間的非彈性碰撞和摩擦消耗系統的能量，從而達到減振的目的[2]。

由于顆粒阻尼復雜的非線性阻尼耗能機理，試驗和仿真成為研究顆粒阻尼的主要手段。Michon等[3]通過實驗探究軟中空顆粒在蜂窩結構中的減振效果，并提出簡化的理論計算模型；Yang等[4]利用功率測量的方法對顆粒阻尼器的耗能特性進行了試驗研究，研究了各因素對阻尼器損耗功率和動態質量的影響；夏兆旺等[5]以填充顆粒的懸臂梁為研究對象，研究了結構阻尼隨顆粒各參數呈非線性變化的特性以及顆粒材料、填充率和金屬盒的位置對系統阻尼的影響；而周宏偉[6]和段勇等[7]在直升機旋翼上使用顆粒阻尼器，試驗表明顆粒阻尼在旋翼非旋轉及低轉速下可以明顯提升系統阻尼，降低旋翼響應水平；Wong等[8]利用PFC3D軟件仿真了顆粒阻尼器損耗功率與振動幅值和振動頻率間的關系，并利用試驗進行了驗證；Xiao等[9]基于FEM-DEM的耦合仿真，對安裝顆粒阻尼器的旋轉齒輪進行了仿真研究，得出了不同尺寸的顆粒的減振效果，并通過實驗進行驗證。

目前，對顆粒阻尼的研究主要集中在對結構的各階模態阻尼的貢獻和減振效果上，對耗散理論及阻尼器內顆粒系統本身的狀態的研究還不夠深入。顆粒物質的特殊性就在于它既可以表現出類似于固態、液態、氣態的特性，又具有豐富而獨特的行為。在一定條件下，它會像液體和氣體一樣流動，表現出一系列流變學行為。但僅有少量學者對顆粒系統不同運動狀態的能量耗散特性進行了研究[10]。并且，對于顆粒物質材料參數對顆粒阻尼器阻尼特性的影響還未被系統的探究。這些在一定程度上阻礙了顆粒阻尼的發展和應用。

近年來，隨著對離散單元法的不斷改進和計算機技術的發展，離散單元法逐漸成為研究顆粒系統的一種重要手段。離散單元法能夠準確得到NOPD中單個顆粒的受力與運動，為分析NOPD中顆粒系統的能量耗散規律提供足夠的信息[11]。

鑒于此，本文采用離散元法對顆粒阻尼器進行了仿真研究，結合阻尼器中顆粒系統的運動狀態探究不同振動參數下顆粒阻尼器的阻尼效果，從細觀尺度上闡釋顆粒系統的能量耗散特性。并且針對影響顆粒阻尼減振特性的內在因素，通過正交試驗，考察了顆粒的材料參數對損耗因子的影響，為顆粒阻尼器的優化設計和應用提供理論依據。

1 仿真模型及參數

1.1 離散單元法

本文使用的顆粒模型是軟球模型，顆粒的法向力被簡化成一個彈簧阻尼元件，切向力被簡化成一個彈簧阻尼元件和一個滑動摩擦元件，并引入彈性系數和阻尼系數等參數來量化彈簧、阻尼器、滑動器的作用，不考慮顆粒表面變形，依據顆粒間法向重疊量和切向位移計算接觸力。所采用的接觸模型是無黏球形顆粒的 Hertz-Mindlin 模型[12]，該模型采用 Hertz 理論計算法向力，采用 Mindlin 方法[13]來計算切向力。設半徑為R1、R2的兩球形顆粒發生彈性接觸，則法向重疊量δn為：

δn=R1+R2-|c1-c2|>0

(1)

式中:R1和R2是兩球的半徑;c1和c2是兩球球心的位置矢量。

顆粒間的接觸面為圓形，則接觸面的半徑a為：

(2)

法向力表示為

(3)

式中:nc為顆粒i球心到顆粒j球心的單位矢量，n=Ri/|Ri|；v為顆粒i球心到顆粒j間的相對運動速度矢量，vnij=vi-vj；kn為法向剛度，根據Hertz接觸理論確定：

(4)

式中:E*和R*分別為有效彈性模量和有效顆粒半徑，可由下式求得：

(5)

(6)

式中:E1、μ1和E2、μ2分別為顆粒1和顆粒2的彈性模量和泊松比。

顆粒間的切向力可表示為：

Fct=-ktδt+dt(vtij·nc)nc

(7)

式中:δt為切向重疊量，即接觸點的切向位移；vtij為接觸點的滑移速度：

vtij=vj-vi+ωj×Rj-ωi×Ri

(8)

kt為切向剛度，由 Mindlin等[13]接觸理論確定，其表達式為：

(9)

G*是有效剪切模量，其值可由下式給定：

(10)

式中:G1、G2是兩顆粒的剪切模量，其值是由彈性模量和泊松比確定的：

(11)

(12)

而對于接觸當中的能量耗散現象，除了庫倫摩擦外，則使用Tsuji等[14]提出的以下非線性黏性阻尼來計算阻尼系數：

(13)

式中:m*為顆粒等效質量；ζ為等效黏性阻尼比，分別由以下兩式求出：

(14)

(15)

式中:mi、mj分別為兩顆粒的質量；e為恢復系數。

1.2 仿真模型

本文建立的3D NOPD 顆粒系統仿真模型如圖1所示，顆粒系統由一個密閉的有機玻璃圓柱容器(D×H=26 mm×78 mm)和裝填在其中的500顆鋼球(d=3.5 mm)組成。仿真采用的模型及材料基于文獻[10]中已驗證的模型，相關參數設置如表1所示。容器受到豎直方向的正弦激勵，仿真的兩個主要控制參數為激振頻率f(15～90 Hz)和激振加速度Γ(1～30 g)。

圖1 DEM顆粒阻尼仿真模型 Fig.1 DEM model of the NOPD

仿真參數數值顆粒材料密度ρ/(kg·m-3)7800顆粒材料彈性模量Ep/Pa206E+009顆粒材料泊松比μp0.3顆粒直徑d/m0.0035總顆粒數量500容器材料密度ρc/(kg·m-3)1190容器材料彈性模量Ec/Pa3.3E+009容器材料泊松比μc0.37恢復系數e0.92靜摩擦因數μs0.3滾動摩阻系數μr0.01時間步Δt/s2.07E-007仿真時間t/s2

2 結果與討論

2.1 顆粒系統運動狀態及相圖

顆粒系統在垂直振動的條件下，會呈現出如類固態[10]、流化、蹦床[15]、對流[16]、萊頓弗羅斯特效應[17]等基于振動顆粒物質流變特性的運動狀態。通過仿真，本文把NOPD中顆粒系統的運動形態分為類固態、局部流化、全局流化、蹦床、對流、萊頓弗羅斯特效應及浮力對流七種。圖2為在不同激振條件下，七種顆粒系統運動形態圖及對應速度矢量圖，從中可以獲得每種運動形態的特性。

當振動加速度低于1.2 g時，顆粒系統表現為類固態，此時容器中的顆粒之間、顆粒與容器壁之間幾乎沒有相對運動，整個阻尼器相當于一個附加質量塊。

當振動加速度達到1.2 g，顆粒物質在“固”、“液”兩種狀態之間的轉變出現。顆粒系統中底部顆粒仍舊沒有發生相對運動，但表面顆粒開始流化，系統進入局部流化狀態，開始表現出類似液態的行為特性。圖2(b)可以看到頂部顆粒層中顆粒的輕微相對運動。

隨著激振頻率和振幅的增加，當局部流化的范圍不斷擴大至整個顆粒系統時，容器中的所有顆粒層的顆粒都發生相對運動，整個顆粒系統處于一個整體流化的狀態，即圖2(c)中的全局流動狀態。

蹦床狀態出現在振幅較大的區域。在蹦床系統中，顆粒就像無彈性的跳躍的小球，在容器中進行雜亂無章的碰撞運動。

對流發生時，不僅顆粒之間存在相對運動，不同顆粒層中的顆粒還會進行顆粒位置的交換，整個顆粒系統就像受熱不均的液體。圖 2(e)對應于頻率為65 Hz時顆粒系統速度矢量圖，可以看到顆粒層中浮動著流動的對流卷。

顆粒的萊頓弗羅斯特效應發生在大振動強度大振幅條件下，顆粒分布出現明顯的Z向密度逆反，即一簇具有六角密排結構的密集顆粒簇被下方一些劇烈運動著的稀少顆粒托起，持續的激勵使下方稀疏的顆粒層不斷發生振蕩，同時上方稠密的顆粒層穩定的浮動在幾乎相同的位置。底部運動劇烈的顆粒相當于蒸汽層，上部顆粒相當于懸浮層[17]。

浮力對流可以看作顆粒萊頓弗羅斯特效應與對流兩種狀態的疊加，當萊頓弗羅斯特現象中懸浮的顆粒層中產生對流運動，可以稱為浮力對流狀態。從圖2(g)中的速度矢量圖可以看到，懸浮的顆粒層中也浮動著流動的對流卷。

此外，還可以通過顆粒系統的勢能曲線看出不同流變行為的特征。本文選取了其中有代表性的四種運動形態隨時間變化的勢能曲線，如圖 3。可以看到，在局部流化狀態時，顆粒系統只有表面顆粒流化，因此勢能變化很小且穩定在一個較小值；同樣波動很小的還有浮力對流狀態，萊頓弗羅斯特效應使大部分顆粒穩定懸浮在特定高度，因此勢能也穩定在略大于局部流化狀態的特定數值；對流狀態時顆粒運動明顯增加，因此勢能變化較大；勢能變化最大的是蹦床狀態，顆粒在整個容器內劇烈運動，勢能波動也相當明顯。

圖3 顆粒系統的勢能曲線 Fig.3 Potential energy curves of granular system in NOPD

為了更清晰直觀地表示顆粒系統出現不同運動狀態的條件，本文繪制了NOPD在不同振動條件下的顆粒系統運動狀態相圖。

從圖4中可以直觀地觀察到顆粒系統運動狀態的轉變，值得注意的是全局流動狀態是系統繼續過渡至其他狀態的基礎。相比其他狀態，對流是一個較不穩定的狀態，它與全局流動的狀態界限不像其他狀態轉化界限那么清晰。萊頓弗羅斯特效應和浮力對流均發生在高振幅高頻率下，但相比萊頓弗萊斯特效應，浮力對流狀態發生在更大的振動強度條件下，即顆粒系統達到產生萊頓弗萊斯特現象的條件后，繼續增大振動強度可使系統轉化為浮力對流狀態。

2.2 顆粒系統損耗因子

損耗因子是衡量系統的阻尼特性并決定其振動能量耗散能力的重要參數。將一個阻尼器的阻尼容量定義為一個完整周期運動中的能量耗散，即

ΔU=Edissipated=fddx

(16)

該表達式通過位移-力平面上的遲滯環路面積得到。如果系統最初(總)的能量表示為Umaximum，那么比阻尼容量D通過下面的比值得到：

(17)

損耗因子η等于在一個周期內每弧度的比阻尼容量。因此

(18)

式中:Edissipated為一個周期內的能量耗散，Emaximum為周期內初始最大能量，近似等于周期內最大動能。在能量分析中，很多重要參數如阻尼比、阻尼系數都與損耗因子有聯系，因此本文引入損耗因子[8, 18]來衡量不同條件下NOPD顆粒系統的阻尼效果。

從圖4中選取四個經歷過較多狀態轉變的頻率，分別繪制出仿真獲得的不同激振頻率下的顆粒系統能量耗散、最大動能及損耗因子曲線。圖5、6可以看到一個周期內的能量耗散隨著振動強度的增加而增加；且頻率越大，增加的幅度越小，大頻率時能量耗散幾乎穩定在一個較小值。周期內最大動能有相似的變化趨勢，但在大頻率大振動強度條件下， 最大動能的值會略微減小。

圖5 不同激振頻率下顆粒系統能量耗散圖6 不同激振頻率下顆粒系統最大動能圖7 不同激振頻率下顆粒系統損耗因子Fig.5TheenergydissipationinacycleversusexcitationintensityunderdifferentfrequencyFig.6themaximumkineticenergyinacycleversusexcitationintensityunderdifferentfrequencyFig.7Lossfactorversusexcitationintensityunderdifferentfrequency

圖7為不同激振頻率下顆粒系統損耗因子隨振動強度變化曲線，對比不同曲線可以發現，不同激振頻率下的顆粒阻尼器表現出最好的阻尼效果時對應的振動強度是不同的。 當系統的激振頻率較小(f=30 Hz或50 Hz)時，損耗因子的最大值出現在振動強度較小的條件下；而當系統的激振頻率較大(f=70 Hz或90 Hz)時，隨著振動強度的變化，損耗因子在經歷波動后會呈不斷上升的趨勢并達到另一個高峰，且損耗因子的最大值出現在振動強度較大的條件下。

結合相圖，可以用顆粒系統運動狀態的轉變來解釋損耗因子數值產生的變化。

當振動強度較小時，顆粒系統表現為類固態。此時NOPD相當于附加質量塊，損耗因子很小。增大振動強度或增大振幅都會使顆粒系統中尤其是位于顆粒系統上層的顆粒運動加劇，使系統狀態由類固態轉化為局部流化，此時顆粒之間，顆粒與容器壁之間的碰撞次數增加，系統耗能增加，但系統最大動能仍然很小，因此損耗因子的數值會有大幅度的提升。

低頻條件下，振動強度的增大意味著振幅的增大，顆粒系統由流化態轉變為蹦床狀態。蹦床狀態中劇烈的碰撞和較多的碰撞機會耗散了很大的能量，而隨著振動條件的變化，頻率的增加和振幅的減小都會使碰撞的劇烈程度顯著減小，使損耗因子的數值減小。同時，損耗因子的數值不會隨著振動強度的增加無限增大，因為此時大部分碰撞是顆粒與容器壁之間的碰撞，這些過于混亂的顆粒運動和碰撞產生倍周期分岔現象，反而不利于系統減振[19]，因此當頻率為30 Hz時，損耗因子的數值很小，幾乎穩定在0.12附近。

相比其他狀態，對流是一個較不穩定的狀態，它與全局流動狀態的界限不像其他狀態轉化界限那樣清晰，出現的范圍也存在一定周期性，隨著振動強度的增大，對流在更高的頻率下出現。因此兩者的損耗因子數值差距也不大，但發生對流運動時，顆粒系統周期變化的方向性給顆粒帶來更多的碰撞機會，這意味著對流狀態時的損耗因子比全局流動狀態略微較大。

萊頓弗羅斯特效應和浮力對流狀態這兩種狀態的發生都是基于顆粒的萊頓弗萊斯特現象，可以說浮力對流是萊頓弗萊斯特現象中的特殊情況，因此都需要滿足高振幅高頻率的振動條件。發生顆粒的萊頓弗萊斯特現象時，顆粒系統上部懸浮著的密集顆粒層被系統底部稀疏的劇烈運動顆粒托起，相比其他狀態，除了顆粒與顆粒之間、顆粒與容器之間碰撞耗散的能量之外，萊頓弗萊斯特現象的發生還需要耗散額外的能量轉化為上部懸浮顆粒層的勢能。因此，這兩種狀態下的損耗因子數值相對較大。但相比萊頓弗萊斯特效應，浮力對流狀態發生在更大的振動強度條件下。這是因為發生萊頓弗羅斯特效應和浮力對流兩種狀態時，顆粒系統中懸浮顆粒層的細觀組成結構是不同的，萊頓弗萊斯特效應中的懸浮層顆粒相當于類固態為六角密排結構，而浮力對流狀態下懸浮層中的顆粒發生的是對流運動，此時顆粒系統需要更多的能量來使懸浮層中的顆粒發生對流運動。類比于類固態和對流狀態時損耗因子數值的關系，浮力對流狀態時懸浮層中顆粒的碰撞與摩擦耗散的能量多于萊頓弗萊斯特效應。因此，當顆粒系統處于浮力對流狀態時，NOPD表現出最好的阻尼效果。

2.3 顆粒材料影響參數正交試驗

NOPD中顆粒材料參數是影響阻尼效果的重要內在原因之一，考慮到外界激勵對于阻尼性態影響重大，為了探索密度、剪切模量、恢復系數、靜摩擦因數、滾動摩擦因數五種顆粒的材料參數對阻尼效果的影響，并綜合考慮計算資源本文在顆粒系統處于七種運動狀態，分別采用正交試驗法進行研究。不僅可以確定各個因素的影響主次，同時對計算資源的需求較低。不同運動狀態對應的外界激勵水平如表2。

表2 不同運動狀態對應的外界激勵水平表Tab.2 External excitation levels corresponding todifferent motion modes table

現以全局流化狀態下正交試驗過程為例，因素和水平的選取如表3所示，選用正交表L16(45)如表4，確定試驗方案。

表3 影響NOPD阻尼效果的因素水平表Tab.3 Factors and their levels used in the orthogonal table

正交試驗方案如表4所示，分別對16組工況進行仿真，并計算顆粒系統的損耗因子，結果如表5所示。其中極差按照下式計算：

R=max{K1,K2,K3,K4}-min{K1,K2,K3,K4}

式中：Ki為任一列水平號為i所對應試驗結果之和，R為極差。

表4 正交表Tab.4 The orthogonal table

由表5分析可知，五個因素在顆粒系統阻尼效果的影響程度上主次是：靜摩擦因數>顆粒密度>剪切模量>滾動摩擦因數>恢復系數。

表5 仿真顆粒系統損耗因子結果表Tab.5 The orthogonal test data and analysis results of theloss factor in the simulation

表6為七種狀態下分別進行正交試驗后得到的極差及影響主次排序。

表6 不同運動狀態下正交試驗結果表Tab.6 The results of orthogonal test data in differentmotion modes

從表6中可以看到， 當顆粒系統處于類固態時，顆粒進入一種密集塊狀的狀態，顆粒的材料參數幾乎對顆粒阻尼器的阻尼性能沒有影響，此結論與前人的研究結論相似[20]。另一個較為特殊的狀態是蹦床狀態，這一狀態下材料參數對顆粒系統阻尼效果有一定影響，但影響效果差距很小。在其他狀態下，材料密度對阻尼效果的影響都比較明顯，這符合前人研究得出的顆粒阻尼減振特性與顆粒材料密度關系的結論[21]。剪切模量和恢復系數對阻尼效果的影響較小。靜摩擦因數和滾動摩擦因數對處于萊頓弗羅斯特效應和浮力對流兩種狀態下的顆粒系統影響較小，但對其余狀態下顆粒系統的阻尼效果有一定影響，原因可能是這兩種狀態下顆粒碰撞及摩擦耗能占總耗能的比例較小。

3 結 論

本文基于離散單元法仿真，結合NOPD中顆粒系統的運動狀態對阻尼效果進行分析，探究不同振動條件對阻尼性能的影響，并通過正交試驗，考察顆粒的材料參數對損耗因子的影響，得到的主要結論如下：

(1)不同激振頻率下，顆粒阻尼器表現出最好的阻尼效果時對應的振動強度是不同的。當系統處于低頻或中頻時，損耗因子的最大值出現在振動強度較小的條件下；而系統處于高頻時，損耗因子的最大值出現在振動強度較大的條件下。

(2)顆粒系統表現出不同運動狀態時的流變學行為及細觀結構導致了阻尼效果的變化，顆粒系統在浮力對流狀態時表現出最優的阻尼效果。

(3)當顆粒系統處于不同的運動狀態時，顆粒的材料參數對系統阻尼的影響程度是不同的。當顆粒系統處于類固態時，顆粒的材料參數幾乎對顆粒阻尼器的阻尼性能沒有影響；蹦床狀態下材料參數對顆粒系統阻尼效果有一定影響，但影響效果差距很小。在其他狀態下，材料密度對阻尼效果的影響都比較明顯，而剪切模量和恢復系數對阻尼效果的影響較小；靜摩擦因數和滾動摩擦因數對處于萊頓弗羅斯特效應和浮力對流兩種狀態下的顆粒系統影響較小。

[ 1 ] 魯正， 呂西林， 閆維明. 顆粒阻尼技術的研究綜述[J]. 振動與沖擊， 2013, 32(7): 1-7.

LU Zheng, Lü Xilin, YAN Weiming. A survey of particle damping techology[J]. Journal of Vibration and Shock, 2013, 32(7): 1-7.

[ 2 ] BAJKOWSKI J M, DYNIEWICZ B, BAJER C I. Damping properties of a beam with vacuum-packed granular damper[J]. Journal of Sound and Vibration, 2015, 341: 74-85.

[ 3 ] MICHON G, ALMAJID A, ARIDON G. Soft hollow particle damping identification in honeycomb structures[J]. Journal of Sound and Vibration, 2013, 332(3): 536-544.

[ 4 ] YANG M Y. Development of master design curves for particle impact dampers[D]. The Pennsylvania State University, 2003.

[ 5 ] 夏兆旺， 單穎春， 劉獻棟. 基于懸臂梁的顆粒阻尼實驗[J]. 航空動力學報， 2007, 22(10): 1737-1741.

XIA Zhaowang, SHAN Yingchun, LIU Xiandong. Experimental research on particle damping of cantilever beam[J]. Journal of Aerospace Power, 2007, 22(10): 1737-1741.

[ 6 ] 周宏偉. 顆粒阻尼及其控制的研究與應用[D]. 南京: 南京航天航空大學, 2008.

[ 7 ] 段勇， 陳前， 林莎. 顆粒阻尼對直升機旋翼槳葉減振效果的試驗[J]. 航空學報, 2009, 30(11): 2113-2118.

DUAN Yong, CHEN Qian, LIN Sha. Experiments of vibration reduction effect of particle damping on helicopter rotor blade[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(11): 2113-2118.

[ 8 ] WONG C, DANIEL M, RONGONG J. Energy dissipation prediction of particle dampers[J]. Journal of Sound and Vibration, 2009, 319(1): 91-118.

[ 9 ] XIAO W, LI J, WANG S, et al. Study on vibration suppression based on particle damping in centrifugal field of gear transmission[J]. Journal of Sound and Vibration, 2016, 366: 62-80.

[10] ZHANG K, CHEN T, WANG X, et al. Rheology behavior and optimal damping effect of granular particles in a non-obstructive particle damper[J]. Journal of Sound and Vibration, 2016, 364: 30-43.

[11] MAO K., WANG M Y, XU Z, et al. DEM simulation of particle damping[J]. Powder Technology, 2004, 142(2): 154-165.

[12] MINDLIN R D. Compliance of elastic bodies in contact[J]. Journal of Applied Mechanics, 1949(16): 259-268.

[13] MINDLIN R D, DERESIEWICZ H. Elastic spheres in contact under varying oblique forces[J]. Journal of Applied Mechanics, 1953(20): 327-344.

[14] TSUJI Y, TANAKA T, ISHIDA T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[J]. Powder Technology, 1992, 71(3): 239-250.

[15] JIANG Z H, WANG Y Y, WU J. Subharmonic motion of granular particles under vertical vibrations[J]. Epl, 2006, 74(3): 417-423.

[16] EHRICHS E, JAEGER H, KARCZMAR G S, et al. Granular convection observed by magnetic resonance imaging[J]. Science, 1995(267): 1632.

[17] ESHUIS P, WEELE K V D, MEER D V D, et al. Granular leidenfrost effect: Experiment and theory of floating particle clusters[J]. Physical Review Letters, 2005( 95): 258001.

[18] YANG B, YICHUNG M. Development of master design curves for particle impact dampers[D]. Doctoral Thesis, The Pennsylvania State University, 2003.

[19] MUJICA N, MELO F. Experimental study of solid-liquid-type transitions in vibrated granular layers and the relation with surface waves[J]. Physical Review E, 2000(63): 011303.

[20] SANCHEZ M, ROSENTHAL G, PUGNALONI L A. Universal response of optimal granular damping devices[J]. Journal of Sound and Vibration, 2012, 331 (20): 4389-4394.

[21] 胡溧. 顆粒阻尼的機理與特性研究[D]. 武漢： 華中科技大學, 2008.