夏兆旺，溫華兵，劉獻(xiàn)棟

(1.江蘇科技大學(xué) 能源與動(dòng)力工程學(xué)院，江蘇 鎮(zhèn)江 212003；2.北京航空航天大學(xué) 交通科學(xué)與工程學(xué)院，北京 100191)

顆粒阻尼是一種被動(dòng)形式阻尼，是在結(jié)構(gòu)內(nèi)部適當(dāng)位置加工出一系列孔腔，也可以是附加在振動(dòng)結(jié)構(gòu)上的獨(dú)立腔體，在這些腔體中填入一定數(shù)量的顆?；蛘叻勰┎牧希w粒常為碳化鎢、陶瓷、鐵、鉛等金屬材料，其形狀可以是球形、橢球形或者不規(guī)則形狀[1-2]。當(dāng)結(jié)構(gòu)發(fā)生振動(dòng)時(shí)，通過(guò)結(jié)構(gòu)孔壁與顆粒的接觸、碰撞將其部分動(dòng)能傳遞給顆粒，顆粒在振動(dòng)結(jié)構(gòu)的孔腔中運(yùn)動(dòng)，進(jìn)而顆粒間、顆粒與結(jié)構(gòu)孔壁間發(fā)生摩擦和碰撞，這些摩擦與碰撞作用將導(dǎo)致能量損耗，振動(dòng)結(jié)構(gòu)的振動(dòng)能量降低，振動(dòng)幅值減小，從而達(dá)到增大系統(tǒng)阻尼的目的[3]。顆粒阻尼技術(shù)具有強(qiáng)非線性，其數(shù)值仿真技術(shù)具有挑戰(zhàn)性。

顆粒阻尼技術(shù)已被廣泛研究，其起源于沖擊阻尼技術(shù)[4]。Panossian等[5-6]對(duì)微顆粒阻尼進(jìn)行了大量研究，研究結(jié)果表明：顆粒阻尼減振效果顯著，在不同種類的顆粒阻尼技術(shù)中顆粒密度越大，效果越好。之后他們將顆粒阻尼用于鋁質(zhì)梁減振以及高頻大振幅的航天飛機(jī)主發(fā)動(dòng)機(jī)進(jìn)口分流葉片上，施加顆粒阻尼后的最大振幅大幅度降低，取得了明顯的減振效果。但對(duì)顆粒阻尼的減振機(jī)理研究不夠深入。

研究表明顆粒阻尼技術(shù)是極端惡劣環(huán)境下結(jié)構(gòu)減振的合適選擇，如飛機(jī)發(fā)動(dòng)機(jī)葉片。近幾十年，研究者廣泛研究了顆粒阻尼技術(shù)在航天主發(fā)動(dòng)機(jī)液氧噴口、航空天線和運(yùn)動(dòng)設(shè)備等領(lǐng)域的減振作用[7-9]。但目前，還未見(jiàn)顆粒阻尼技術(shù)應(yīng)用于旋轉(zhuǎn)結(jié)構(gòu)減振的理論研究報(bào)道。

為了研究顆粒阻尼技術(shù)在葉片結(jié)構(gòu)減振上應(yīng)用的可行性，同時(shí)便于計(jì)算，本文將分析顆粒阻尼技術(shù)在平板葉片結(jié)構(gòu)減振上的應(yīng)用研究。數(shù)值仿真研究能加深對(duì)顆粒阻尼減振復(fù)雜機(jī)理的理解，也能降低顆粒阻尼技術(shù)應(yīng)用的設(shè)計(jì)成本。本文將研究基于離散元和有限元耦合算法的帶顆粒阻尼的旋轉(zhuǎn)平板葉片結(jié)構(gòu)動(dòng)力學(xué)模型，分析不同參數(shù)對(duì)顆粒阻尼減振效果的影響規(guī)律。

1 耦合仿真算法

本文針對(duì)帶顆粒阻尼器的平板葉片結(jié)構(gòu),將擅長(zhǎng)計(jì)算復(fù)雜連續(xù)結(jié)構(gòu)的有限元法與擅長(zhǎng)計(jì)算離散顆粒運(yùn)動(dòng)的離散元法結(jié)合起來(lái)，提出了基于有限元—離散元的耦合仿真算法。應(yīng)用DMAP語(yǔ)言對(duì)工程有限元軟件NASTRAN進(jìn)行二次開(kāi)發(fā)，通過(guò)MATLAB將有限元法與離散元法進(jìn)行耦合。

圖1 仿真算法流程圖

帶顆粒阻尼器旋轉(zhuǎn)結(jié)構(gòu)的振動(dòng)特性計(jì)算流程如圖1所示，主要通過(guò)MATLAB調(diào)用相應(yīng)的有限元和離散元計(jì)算結(jié)果進(jìn)行循環(huán)計(jì)算，其整個(gè)流程分為：①建立帶顆粒阻尼器結(jié)構(gòu)的有限元模型，并根據(jù)顆粒大小、相對(duì)速度及顆粒阻尼器孔腔尺寸，通過(guò)MATLAB隨機(jī)生成各個(gè)顆粒的初始位置，并設(shè)置顆粒初始速度為零。 ②利用有限元軟件NASTRAN計(jì)算結(jié)構(gòu)在激振力作用下的響應(yīng)，提取結(jié)構(gòu)相應(yīng)有限元節(jié)點(diǎn)的位移、速度、加速度等信息。其中，提取結(jié)構(gòu)上合適位置的加速度用于研究結(jié)構(gòu)振動(dòng)特性；提取孔腔周圍各點(diǎn)的位移和速度用于計(jì)算孔腔與顆粒間的接觸力。③通過(guò)MATLAB調(diào)用孔腔周圍各點(diǎn)的位移和速度，根據(jù)位移關(guān)系判斷顆粒與孔腔間是否接觸。若接觸，即根據(jù)疊合量和相對(duì)速度求出顆粒間的接觸力，計(jì)算一個(gè)步長(zhǎng)內(nèi)各個(gè)顆粒在接觸力作用下的位移、速度，更新各個(gè)顆粒的位移和速度，并計(jì)算結(jié)構(gòu)的孔壁各個(gè)位置受到顆粒的接觸力。④通過(guò)MATLAB將接觸力賦予NASTRAN軟件計(jì)算程序，計(jì)算結(jié)構(gòu)在激振力和接觸力共同作用下一個(gè)步長(zhǎng)的響應(yīng)，并提取相應(yīng)點(diǎn)的位移、速度、加速度等信息。⑤根據(jù)計(jì)算時(shí)間判斷循環(huán)是否結(jié)束。如果結(jié)束,輸出結(jié)構(gòu)相應(yīng)點(diǎn)的加速度，一般需計(jì)算在激振力作用下十個(gè)以上周期響應(yīng)，否則繼續(xù)循環(huán)執(zhí)行第③，④，⑤步驟，直到程序結(jié)束。

1.1 接觸力學(xué)模型

圖2 顆粒接觸模型

離散元法可以模擬運(yùn)動(dòng)在顆粒集合中傳播的過(guò)程。顆粒運(yùn)動(dòng)必然會(huì)導(dǎo)致顆粒之間的相互碰撞，顆粒之間也必然有接觸力產(chǎn)生。球狀顆粒間產(chǎn)生接觸力的條件是兩球心之間的距離小于兩球半徑之和。兩球狀顆粒接觸時(shí)的相對(duì)運(yùn)動(dòng)有三種：在連心線方向上的法向運(yùn)動(dòng)；在兩接觸面的相對(duì)切向運(yùn)動(dòng)；兩球的相對(duì)轉(zhuǎn)動(dòng)。本文在計(jì)算顆粒間的接觸力時(shí)，只考慮法向和切向接觸力，并分別采用不同的接觸力學(xué)模型計(jì)算，因?yàn)轭w粒接觸面積很小，由兩球球心方向相對(duì)轉(zhuǎn)動(dòng)引起的力矩可忽略不計(jì)。顆粒單元之間的接觸可以簡(jiǎn)化為彈簧-阻尼器-滑塊模型，如圖2所示。彈簧代表單元的彈性，阻尼器代表單元的粘性，用滑塊來(lái)表示單元之間的摩擦。

赫茲公式只適用于無(wú)摩擦的彈性接觸，所以，它只能對(duì)顆粒接觸的法向應(yīng)力進(jìn)行計(jì)算。赫茲黏彈性模型表示的法向作用力為

(1)

式中:kn為法向剛度系數(shù)，cn為法向阻尼系數(shù)，δn為疊合量。

如果法向位移恒定，法向力和接觸面半徑都將保持不變，則切向力隨實(shí)際的切向位移而變。但法向位移與切向位移往往同時(shí)發(fā)生變化，此時(shí)可以通過(guò)切向的位移增量與接觸力增量間的關(guān)系得到切向力，其粘彈性力學(xué)模型為

Fs(t)=Fns(t-Δt)+ksvsΔt+csvs

(2)

切向力還需滿足庫(kù)侖定律，即:

(3)

式中:ks為切向剛度系數(shù)，vs為切向速度，cs為切向阻尼系數(shù)。Fs(t)為t時(shí)刻接觸兩體間的切向作用力；Fns(t-Δt)為t-Δt時(shí)刻接觸兩體間的切向作用力；Δt為計(jì)算時(shí)間步長(zhǎng)；ks為接觸兩體的切向剛度系數(shù)，可由經(jīng)驗(yàn)或?qū)嶒?yàn)確定，也可由下式計(jì)算

帶顆粒阻尼器結(jié)構(gòu)在旋轉(zhuǎn)運(yùn)動(dòng)時(shí)，顆粒除了受到顆粒間及顆粒與孔腔間的接觸力外，還受到離心力作用。其離心力的大小為

P=mω2r

(4)

其中:m為顆粒的質(zhì)量；ω為結(jié)構(gòu)的轉(zhuǎn)速；r為顆粒的離心半徑。

1.2 顆粒接觸時(shí)步

在計(jì)算中，時(shí)步必須取得很小，以保證在一個(gè)時(shí)步內(nèi)每個(gè)顆粒的動(dòng)量只能傳播給其臨近的顆粒，以及顆粒加速度是近似恒定的。實(shí)際中，通常根據(jù)瑞利波的傳播特性確定時(shí)步。根據(jù)瑞利波沿著顆粒表面的傳播速度得到臨界時(shí)間步長(zhǎng)為

(5)

式中:G為剪切模量，ρ為顆粒密度，γ為泊松比，Rmin為顆粒最小半徑。

根據(jù)瑞利波確定的時(shí)步一般能滿足數(shù)值仿真收斂的要求。對(duì)于顆粒之間相對(duì)速度很大的情況，當(dāng)與一個(gè)顆粒接觸的顆粒數(shù)大于等于4時(shí)，此時(shí)這個(gè)顆粒的接觸時(shí)間取公式中的20%較為合適；當(dāng)與一個(gè)顆粒接觸的顆粒數(shù)小于4時(shí)，此時(shí)這個(gè)顆粒的接觸時(shí)間取公式中的40%較為合適。事實(shí)上，盡管很多學(xué)者對(duì)合理時(shí)步的選取做了大量工作，但由于時(shí)步的選取與很多因素有關(guān)，因此在實(shí)際計(jì)算時(shí)還應(yīng)結(jié)合實(shí)際情況通過(guò)試算來(lái)確定。

1.3 顆粒運(yùn)動(dòng)方程

離散元法的原理和計(jì)算過(guò)程是根據(jù)單元間力的相互作用和牛頓運(yùn)動(dòng)定律，采用動(dòng)態(tài)松弛法進(jìn)行循環(huán)迭代計(jì)算，在每一個(gè)時(shí)步都更新單元的位置，并遍及整個(gè)單元集合。通過(guò)對(duì)系統(tǒng)中每個(gè)單元的微觀運(yùn)動(dòng)過(guò)程進(jìn)行跟蹤研究，考慮每個(gè)單元對(duì)結(jié)構(gòu)系統(tǒng)的作用力，即可得到整個(gè)系統(tǒng)的宏觀運(yùn)動(dòng)規(guī)律。對(duì)于球形單元在任意坐標(biāo)軸、任意時(shí)刻t，考慮每一單元受力作用后產(chǎn)生的平動(dòng)和轉(zhuǎn)動(dòng)。由牛頓第二定律得到

(6)

式中:xi(t)為平動(dòng)位移、∑F(t)為合力,mi為質(zhì)量，∑M(t)為合力矩，θi(t)為角位移，Ii為轉(zhuǎn)動(dòng)慣量。

由中心差分法可得t+Δt/2時(shí)刻顆粒的速度和加速度

(7)

(8)

這樣，經(jīng)過(guò)Δt后單元i就移動(dòng)到一個(gè)新的位置，并產(chǎn)生新的接觸力和接觸力矩，計(jì)算其所受的合力∑F(t+Δt/2)和合力矩∑M(t+Δt/2)后返回式(8)計(jì)算，重復(fù)這一過(guò)程，即可得到每個(gè)單元以及整個(gè)顆粒群體的運(yùn)動(dòng)特性。

2 旋轉(zhuǎn)平板葉片結(jié)構(gòu)減振仿真結(jié)果

平板葉片結(jié)構(gòu)選擇密度為78 000 kg/m3的鑄鋼，其尺寸為250 mm×60 mm×4 mm。顆粒之間及顆粒與平板葉片結(jié)構(gòu)之間的彈性恢復(fù)系數(shù)都為0.78，顆粒選擇直徑為2 mm的鋼球和20目碳化鎢兩種材料。

通過(guò)模態(tài)實(shí)驗(yàn)確定懸臂平板葉片結(jié)構(gòu)的打孔位置。平板葉片結(jié)構(gòu)的前三階模態(tài)振型如圖3所示，模態(tài)振型顯示了平板葉片結(jié)構(gòu)的最大和最小位移的位置。為研究顆粒阻尼器的孔腔布置方式和減振效果之間的關(guān)系，顆粒阻尼器的孔腔布置方式選擇如圖4所示。

圖3 懸臂平板葉片結(jié)構(gòu)模態(tài)圖

平板葉片結(jié)構(gòu)布置3個(gè)孔腔，正弦激振力的幅值為0.8 N，激振力的頻率范圍分別選擇為30～60 Hz/300～365 Hz/460～550 Hz以覆蓋前三階模態(tài)頻率，平板葉片的旋轉(zhuǎn)速度分別選擇0 r/min和1 000 r/min，激振力位置如圖5所示。顆粒填充率分別為0%(無(wú)顆粒)和70%。圖5～7顯示了旋轉(zhuǎn)平板葉片自由端的響應(yīng)均方根值與激振力頻率之間的關(guān)系及不同填充顆粒對(duì)減振效果的影響。結(jié)果表明：本文提出的方法用于計(jì)算帶顆粒阻尼的旋轉(zhuǎn)板結(jié)構(gòu)響應(yīng)是可行的。

圖4 旋轉(zhuǎn)平板葉片結(jié)構(gòu)的孔腔布置方式

圖5 旋轉(zhuǎn)速度對(duì)平板葉片結(jié)構(gòu)一階響應(yīng)的影響

圖6 旋轉(zhuǎn)速度對(duì)平板葉片結(jié)構(gòu)二階響應(yīng)的影響

圖7 旋轉(zhuǎn)速度對(duì)平板葉片結(jié)構(gòu)三階響應(yīng)的影響

圖8 豎向各組打孔位置示意圖

下面仿真分析平板葉片在旋轉(zhuǎn)狀態(tài)下帶顆粒阻尼器的減振效果。為研究孔洞位置對(duì)振動(dòng)特性的影響，孔洞位置分豎向打孔和橫向打孔兩種情況，豎向打孔又分三組位置——A組、B組、C組,如圖8所示。豎向A組的位置：中間孔腔軸線與固定端的距離為220 mm；豎向B、C組的位置分別是：中間孔腔底部與固定端的距離分別為120 mm和60 mm。無(wú)論橫向還是豎向打孔，每組都是打3個(gè)孔。各種布置的顆粒阻尼器填充率均為70%，填充顆粒為2 mm鋼球，正弦激振力幅值為0.8 N。旋轉(zhuǎn)平板葉片的固有頻率隨轉(zhuǎn)速改變而改變，本文激振力的頻率分別為各轉(zhuǎn)速下的第一階、第二階、第三階固有頻率，旋轉(zhuǎn)葉片在激振力作用下的響應(yīng)曲線分別如圖9～圖11所示。

從圖9～圖11可以看出：平板葉片一、二階響應(yīng)在轉(zhuǎn)速較低時(shí)較小，但隨著轉(zhuǎn)速的增加而增加，當(dāng)轉(zhuǎn)速達(dá)到一定值時(shí),轉(zhuǎn)速對(duì)響應(yīng)的影響趨于穩(wěn)定，其原因是：在轉(zhuǎn)速低時(shí)顆粒受到的離心力小，顆粒發(fā)生積聚后彼此之間壓力還小、約束小，顆粒從孔壁獲得動(dòng)能尚能在聚積的顆粒堆中運(yùn)動(dòng)、并與其它顆粒碰撞、摩擦，使得更多顆粒運(yùn)動(dòng)起來(lái)通過(guò)碰撞、摩擦消耗掉動(dòng)能。然而隨著轉(zhuǎn)速的提高，由于離心力變大使得顆粒之間的壓緊力增加，顆粒之間的約束也增加，與孔壁直接接觸的顆粒獲得的動(dòng)能減少，顆粒之間的碰撞和動(dòng)能的傳遞減少，相對(duì)運(yùn)動(dòng)的顆粒減少，因此在耗能機(jī)制中碰撞耗能大幅減少，同時(shí)相對(duì)運(yùn)動(dòng)的顆粒數(shù)量減少，摩擦消耗也減少，故此表現(xiàn)出的阻尼減少，振動(dòng)響應(yīng)變大；在轉(zhuǎn)速達(dá)到一定值時(shí)，顆粒在較大離心力作用下聚積在一起，顆粒的約束大，顆粒之間幾乎無(wú)碰撞，顆粒與旋轉(zhuǎn)平板葉片的相對(duì)速度接近零，只能依靠結(jié)構(gòu)在外界激勵(lì)下發(fā)生彈性變形處顆粒間的相對(duì)運(yùn)動(dòng)和摩擦消耗能量，此時(shí)消耗的能量受轉(zhuǎn)速影響很小，由于顆粒之間摩擦耗能大于結(jié)構(gòu)內(nèi)耗能，因此仍表現(xiàn)出一定的阻尼效果。此外，還可從圖中看出，平板葉片的第三階響應(yīng)曲線隨轉(zhuǎn)速的增加先減小后增加，然后趨于平穩(wěn)，其原因是葉片在三階共振時(shí)系統(tǒng)的能量主要依靠顆粒及顆粒與孔壁間的摩擦消耗掉，隨著轉(zhuǎn)速由零逐漸增大，離心力使顆粒間的相互作用力增加，顆粒之間的摩擦作用更加明顯，但隨著轉(zhuǎn)速的進(jìn)一步加大，巨大的離心力使得顆粒間的相互運(yùn)動(dòng)很小，系統(tǒng)的能量消耗也隨之降低。

圖9 平板葉片旋轉(zhuǎn)狀態(tài)下響應(yīng)與轉(zhuǎn)速的關(guān)系(一階)

圖10 平板葉片旋轉(zhuǎn)狀態(tài)下響應(yīng)與轉(zhuǎn)速的關(guān)系(二階)

圖11 平板葉片旋轉(zhuǎn)狀態(tài)下響應(yīng)與轉(zhuǎn)速的關(guān)系(三階)

從圖9～圖11還可看出，各種不同形式的顆粒阻尼器布置對(duì)各階振動(dòng)的抑制規(guī)律基本一致，即在轉(zhuǎn)速較低時(shí)抑制效果好，隨著轉(zhuǎn)速提高抑制效果迅速減弱，在轉(zhuǎn)速達(dá)到一定值時(shí)，抑制效果趨于穩(wěn)定，基本不隨轉(zhuǎn)速提高而改變，但在各種轉(zhuǎn)速下抑制效果有所不同，總體來(lái)說(shuō)橫向布置顆粒阻尼器的減振效果相對(duì)較好。

3 結(jié) 論

本文針對(duì)帶顆粒阻尼器的旋轉(zhuǎn)平板葉片結(jié)構(gòu)提出了基于離散元和有限元的耦合仿真算法。仿真結(jié)果表明：①所提出的耦合仿真算法在很寬的頻帶范圍內(nèi)都是可行的；②帶顆粒阻尼器的旋轉(zhuǎn)平板葉片結(jié)構(gòu)的減振效果受平板葉片的轉(zhuǎn)速影響很大，帶顆粒阻尼器的旋轉(zhuǎn)平板葉片結(jié)構(gòu)的第一、二階振動(dòng)響應(yīng)隨轉(zhuǎn)速的提高而增加，即阻尼作用有所降低。第三階振動(dòng)響應(yīng)先降低，在轉(zhuǎn)速為600 r/min時(shí)達(dá)到最低，然后振動(dòng)響應(yīng)隨轉(zhuǎn)速提高迅速增加，在轉(zhuǎn)速達(dá)到1 000 r/min后趨于平穩(wěn)。

[1]夏兆旺,溫華兵,劉獻(xiàn)棟.顆粒阻尼器結(jié)構(gòu)振動(dòng)特性耦合算法仿真與試驗(yàn)[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),2011,42(8):26-29.

XIA Zhao-wang,WEN Hua-bing，LIU Xian-dong.Coupling algorithm simulation and experiment of structure vibration characteristics for particle dampers[J].Transactions of the Chinese Society for Agricultural Machinery,2011,42(8):26-29.

[2]Xia Z W，Liu X D，Shan Y C.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-204.

[3]Zheng L，Xilin L，Masri S F.Studies of the performance of particle dampers under dynamic loads [J].Journal of Sound and Vibration，2010,329:5415-5433

[4]Simonian S S.Particle beam damper[C].Proceedings of SPIE Conference on Passive Damping，SPIE，1995:149-160.

[5]Panossian H V，Johnson V J，Rogers L.Structural damping enhancement via non-obstructive particle damping technique [J].Journal of Vibration and Acoustics，1992，114(1):101-105.

[6]Panossian H V.Structural damping enhancement via non-obstructive particle damping technique[J].ASME Journal of Vibration and Acoustics，1992,114(1):101-105.

[7]Wong C X，Rongong J A.Micro-model characterisation and application of particle dampers to vibrating structures [C].47th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Confer，Newport，Rhode Island 2006:211-224.

[8]Mao K，Wang M，Ding H，et al，Simulation and characterization of particle damping in transient vibrations[J].Journal of Vibra-tion and Acoustics，2004,126 (2):202-211

[9]Papalou A，Masri S F.An experimental investigation of particle damper under harmonic excitation[J].Journal of Vibration and Control，1998,4(4):361-379.