馬皓曄，李琳,2，范雨,2,*，吳亞光

1. 北京航空航天大學(xué) 能源與動(dòng)力工程學(xué)院，北京 100083 2. 北京航空航天大學(xué) 航空發(fā)動(dòng)機(jī)結(jié)構(gòu)強(qiáng)度北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083

結(jié)構(gòu)振動(dòng)普遍存在于航空、航天、船舶等工程領(lǐng)域[1]中。過(guò)大的振動(dòng)會(huì)導(dǎo)致結(jié)構(gòu)的高周疲勞甚至失效，而增加阻尼是降低振動(dòng)的有效方式。雖然對(duì)新型阻尼器的研究長(zhǎng)盛不衰(如壓電分支阻尼技術(shù)[2-5])，但基于摩擦耗能原理的干摩擦阻尼器由于具有結(jié)構(gòu)簡(jiǎn)單、可靠性高、減振效果明顯且受溫度影響較小等優(yōu)點(diǎn)，仍是最常用且最受關(guān)注的阻尼技術(shù)之一[6-8]。在現(xiàn)役航空發(fā)動(dòng)機(jī)中，就有多種形式的干摩擦阻尼器，如葉冠、葉片凸肩、緣板阻尼器、阻尼環(huán)等[9-12]。隨著未來(lái)飛行器對(duì)更高的推重比的追求，將進(jìn)一步大量采用薄壁結(jié)構(gòu)，這必然會(huì)導(dǎo)致更為突出的結(jié)構(gòu)振動(dòng)問(wèn)題[9]。然而針對(duì)適用于一般薄壁結(jié)構(gòu)的干摩擦阻尼器的研究尚不多見(jiàn)。

設(shè)計(jì)干摩擦阻尼器的關(guān)鍵技術(shù)之一是發(fā)展含接觸環(huán)節(jié)的非線(xiàn)性動(dòng)力學(xué)模型及其求解方法。如果主要關(guān)注的是一類(lèi)干摩擦系統(tǒng)的普適動(dòng)力學(xué)行為，則常用集總參數(shù)模型，即將結(jié)構(gòu)用少數(shù)幾個(gè)質(zhì)量和彈簧的組合表征。研究人員已經(jīng)用這類(lèi)模型探索了葉根阻尼器[13]、緣板阻尼器[14-16]、針對(duì)齒輪扭振的單個(gè)[17]和多個(gè)[18]阻尼環(huán)。集中參數(shù)模型的自由度很少，數(shù)值問(wèn)題不突出，因此適合用于機(jī)理性研究或算法的合理有效性研究[19-20]。

如果主要關(guān)注的是具體的干摩擦阻尼器及其對(duì)所施加結(jié)構(gòu)的減振效果(尤其是針對(duì)多個(gè)高階模態(tài)時(shí))，則應(yīng)考慮使用高保真有限元模型。例如，張大義等[21]基于實(shí)體有限元模型建立了緣板阻尼結(jié)構(gòu)的設(shè)計(jì)流程和優(yōu)化方法，并給出了結(jié)構(gòu)設(shè)計(jì)合理性的評(píng)估參數(shù)。Petrov和Ewins[22]利用高保真有限元模型分析了帶冠葉片的強(qiáng)迫響應(yīng)，并對(duì)屋頂形阻尼器(Cottage Roof Dampers, CRDs)與分離式阻尼器(Split Dampers, SDs)的減振性能進(jìn)行了比較[23]，并說(shuō)明相比于CRDS，SDs對(duì)質(zhì)量參數(shù)的變化不敏感。Charleux等[24]將榫接葉盤(pán)結(jié)構(gòu)強(qiáng)迫響應(yīng)的試驗(yàn)與仿真結(jié)果進(jìn)行對(duì)比，說(shuō)明了利用高保真模型來(lái)預(yù)測(cè)帶接觸系統(tǒng)的振動(dòng)特性是可行、可信的。高保真模型的自由度數(shù)一般較大，使得直接進(jìn)行時(shí)域積分變得非常耗時(shí)，這就需要發(fā)展減縮模型和高效的頻域求解算法。目前常用動(dòng)力凝縮和子結(jié)構(gòu)方法減縮線(xiàn)性自由度的規(guī)模[25-26]。對(duì)于非線(xiàn)性自由度，可以在諧波平衡的框架內(nèi)利用周期對(duì)稱(chēng)性對(duì)葉盤(pán)等循環(huán)周期結(jié)構(gòu)上的接觸點(diǎn)數(shù)目進(jìn)行減縮[27]；或?qū)τ谝话憬Y(jié)構(gòu)，Krack等[28-29]還發(fā)展了非線(xiàn)性模態(tài)減縮技術(shù)，用以計(jì)算失諧榫接葉盤(pán)的強(qiáng)迫響應(yīng)。目前最常用的頻域求解方法是諧波平衡法，即通過(guò)傅里葉變換將對(duì)非線(xiàn)性常微分方程的求解轉(zhuǎn)換為對(duì)非線(xiàn)性代數(shù)方程組的求解。即使針對(duì)集總參數(shù)模型，其計(jì)算效率遠(yuǎn)高于時(shí)間積分法[30](CPU時(shí)間只是后者的約1%),但隨著自由度數(shù)的增加，諧波平衡法也面臨多個(gè)數(shù)值問(wèn)題，其中的一個(gè)主要問(wèn)題是用數(shù)值差分求雅克比矩陣的誤差導(dǎo)致牛頓-拉夫遜迭代難以收斂[31]。

本文提出一種適用于一般薄壁結(jié)構(gòu)的波紋形干摩擦阻尼器，具有適用性強(qiáng)、正壓力調(diào)節(jié)方便、易于安裝等特點(diǎn)；并利用高保真有限元模型計(jì)算其頻響特性,展示其阻尼性能隨安裝螺栓壓緊力的變化規(guī)律。在計(jì)算方法上，采用了加速動(dòng)態(tài)拉格朗日法對(duì)強(qiáng)迫響應(yīng)進(jìn)行預(yù)測(cè)，并給出了解析形式的雅克比矩陣，使算法的精度和速度得以大幅提升。通過(guò)集中參數(shù)模型的算例說(shuō)明了理論推導(dǎo)及程序編寫(xiě)的合理與正確性。利用一個(gè)四邊固支的平板結(jié)構(gòu)來(lái)說(shuō)明波紋形摩擦片的阻尼性能。

1 波紋形摩擦阻尼片

本文所考慮的波紋形干摩擦阻尼片如圖1所示。阻尼片的2個(gè)接觸面對(duì)稱(chēng)分布在安裝面兩側(cè)，通過(guò)振動(dòng)時(shí)接觸面與底部承載結(jié)構(gòu)之間的相對(duì)運(yùn)動(dòng)產(chǎn)生的摩擦來(lái)達(dá)到減振的目的。安裝面利用螺栓與螺母進(jìn)行固定。

圖1 波紋形干摩擦阻尼器示意圖
Fig.1 Illustration of corrugated frictional patch

接觸面的正壓力是干摩擦阻尼器設(shè)計(jì)時(shí)的關(guān)鍵參數(shù)，對(duì)干摩擦阻尼器的性能有著顯著的影響。當(dāng)正壓力過(guò)小時(shí)，接觸面產(chǎn)生的切向力較小，導(dǎo)致由摩擦損耗的能量較??；當(dāng)正壓力過(guò)大時(shí)，由于接觸面之間的相對(duì)滑移量減小，使得減振性能下降；當(dāng)正壓力增大到使得接觸面達(dá)到“完全阻滯”狀態(tài)時(shí)，此時(shí)干摩擦阻尼器只有微小的調(diào)節(jié)固有頻率的作用，而不具備阻尼作用。對(duì)于圖1(a)所示的阻尼器，可以通過(guò)改變螺栓的擰緊程度來(lái)調(diào)節(jié)接觸面上的正壓力。

在對(duì)波紋形阻尼器的減振特性進(jìn)行分析時(shí)，采用了如圖1(b)所示的有限元模型。其中T1、W1分別代表阻尼片的厚度與寬度。摩擦片的輪廓線(xiàn)(在xOy平面的投影)關(guān)于中心(安裝點(diǎn))對(duì)稱(chēng)，共由3段曲線(xiàn)構(gòu)成：中間段為長(zhǎng)為L(zhǎng)1的直線(xiàn)，左右2段均為余弦曲線(xiàn),滿(mǎn)足如下關(guān)系：

(1)

式中：L2為半波長(zhǎng);H1為波的峰值。直線(xiàn)段分別與兩邊的余弦曲線(xiàn)段相切。工作過(guò)程中會(huì)通過(guò)擰緊安裝面兩側(cè)的螺母來(lái)減小非摩擦界面處的相對(duì)運(yùn)動(dòng)，此時(shí)安裝面近似于固支狀態(tài)，因此在有限元分析中，螺栓與螺母結(jié)構(gòu)可用圖1(b)中位于安裝面區(qū)域內(nèi)連接摩擦片與主體結(jié)構(gòu)的實(shí)體單元來(lái)代替。接觸節(jié)點(diǎn)位于阻尼片下側(cè)的接觸面上，它們與主體結(jié)構(gòu)上對(duì)應(yīng)的節(jié)點(diǎn)重合。

接觸節(jié)點(diǎn)之間的相對(duì)運(yùn)動(dòng)相比于整體結(jié)構(gòu)的振動(dòng)幅值較小，因此可以采用點(diǎn)對(duì)點(diǎn)式的接觸單元來(lái)描述微滑移運(yùn)動(dòng)，如圖2所示。該模型又被稱(chēng)為三維接觸運(yùn)動(dòng)模型，其中FN代表法向正壓力，由接觸對(duì)之間沿法向的相對(duì)位移決定，因此其能夠考慮由于自身形變所導(dǎo)致正壓力改變；FT1與FT2分別代表接觸平面上2個(gè)正交方向的摩擦力，在計(jì)算中認(rèn)為它們是互不影響的[32]。該模型能夠完備地描述三維空間中接觸對(duì)之間的運(yùn)動(dòng)關(guān)系，在對(duì)高保真有限元模型研究時(shí)應(yīng)用較為廣泛[22-23,26]。

圖2 三維接觸運(yùn)動(dòng)模型示意圖
Fig.2 Illustration of 3D contact movement model

2 干摩擦結(jié)構(gòu)響應(yīng)的快速預(yù)測(cè)

為了展示波紋形干摩擦阻尼片的減振效果，用頻響函數(shù)的峰值作為評(píng)價(jià)指標(biāo)。為此，采用了基于時(shí)頻轉(zhuǎn)換(Alternating Frequency-Time domain, AFT)的速度型動(dòng)態(tài)拉格朗日法[26](Dynamic Lagrange Frequency-Time method， DLFT)。相比于以遲滯彈簧模型為基礎(chǔ)的強(qiáng)迫響應(yīng)預(yù)測(cè)算法[27]，DLFT的主要優(yōu)點(diǎn)為：① 不需要確定切法向彈簧剛度的值，同時(shí)計(jì)算結(jié)果對(duì)引入的懲罰系數(shù)的值不敏感[33]；② 時(shí)域內(nèi)不需要通過(guò)迭代來(lái)使得非線(xiàn)性力滿(mǎn)足周期條件[34]。在此基礎(chǔ)上，引入了解析形式的雅克比矩陣使得算法的計(jì)算效率與收斂穩(wěn)定性大幅提升。

2.1 高階諧波平衡法求解非線(xiàn)性動(dòng)力學(xué)方程

對(duì)于帶有干摩擦阻尼的非線(xiàn)性系統(tǒng)，其離散形式的動(dòng)力學(xué)方程可以寫(xiě)為

(2)

式中：M、C和K分別為質(zhì)量、阻尼和剛度矩陣；qex(t)為外激振力向量，在本文中僅考慮其為時(shí)間的簡(jiǎn)諧形式；qnl(t)代表作用于摩擦界面處的非線(xiàn)性干摩擦力，包括了切向與法向分量；u(t)為位移向量；上標(biāo)(·)代表對(duì)時(shí)間的一次求導(dǎo)。直接在時(shí)域內(nèi)求解式(2)這樣的2階非線(xiàn)性微分方程組將會(huì)耗費(fèi)大量時(shí)間，尤其當(dāng)自由度較多且求解的是穩(wěn)態(tài)響應(yīng)時(shí)，計(jì)算時(shí)間成本將迅速增加。故本文在頻域內(nèi)采用高階諧波平衡法求解式(2)的穩(wěn)態(tài)強(qiáng)迫響應(yīng)。

高階諧波平衡法的出發(fā)點(diǎn)是將位移、外激勵(lì)和非線(xiàn)性力表示成傅里葉級(jí)數(shù)，并截取前Nh項(xiàng)。因此，位移項(xiàng)u(t)可寫(xiě)為

(3)

(4)

(5)

(6)

由于ejnωt≠0，因此只能前面括號(hào)里的向量為零，從而獲得一組非線(xiàn)性代數(shù)方程組：

(7)

式中：H(n)為第n階諧波的動(dòng)剛度矩陣，即

H(n)=-(nω)2M+jnωC+K

(8)

在解該非線(xiàn)性代數(shù)方程組時(shí)，一般常用“牛頓-拉夫遜”迭代法以及其改進(jìn)方法[15]。然而隨著方程自由度的提高，計(jì)算效率與穩(wěn)定性都會(huì)明顯下降。注意到相對(duì)于整個(gè)結(jié)構(gòu)系統(tǒng)，接觸面只占很小一部分，因此本文進(jìn)行自由度凝縮，以減少迭代時(shí)方程的自由度。首先，將式(7)中的向量和矩陣分成線(xiàn)性自由度(下標(biāo)l)與非線(xiàn)性自由度(下標(biāo)nl)：

(9)

再通過(guò)凝聚，僅保留非線(xiàn)性自由度：

(10)

(11)

這樣，總自由度從(Nl+Nnl)×(Nh+1)變?yōu)镹nl×(Nh+1)，其中Nl和Nnl分別為線(xiàn)性自由度數(shù)與非線(xiàn)性自由度數(shù)。

(12)

式中：S(n)=[D(n)]-1。對(duì)于任意一組接觸點(diǎn)，引入相對(duì)位移：

(13)

聯(lián)立方程式(12)與式(13)，則動(dòng)力學(xué)方程組可最終表示為相對(duì)位移的形式：

(14)

其中：

(15)

由此，方程組的自由度數(shù)進(jìn)一步降低到原來(lái)的一半，即Nnl×(Nh+1)/2。在對(duì)式(14)進(jìn)行求解時(shí)，一般將其寫(xiě)為收斂殘差的形式，即

(16)

2.2 速度型動(dòng)態(tài)拉格朗日法求解非線(xiàn)性力

(17)

并用如下關(guān)系計(jì)算非線(xiàn)性力：

(18)

式(17)可以看做是在式(14)的基礎(chǔ)上，分別在切向(下標(biāo)T)引入了一個(gè)關(guān)于相對(duì)速度的懲罰項(xiàng)以及在法向(下標(biāo)N)引入了一個(gè)關(guān)于相對(duì)位移的懲罰項(xiàng)，其中ε為懲罰系數(shù)。

(19)

利用傅里葉逆變換，將式(19)的每一項(xiàng)都變換到時(shí)域，此時(shí)每個(gè)自由度的運(yùn)動(dòng)都被離散至Ns個(gè)時(shí)間點(diǎn)上，其中任意第l個(gè)時(shí)間點(diǎn)下的拉格朗日乘子λ(l)可表示為

λ(l)=λu(l)-λx(l)

(20)

注意此時(shí)實(shí)際上已知的是λu，λ和λx暫時(shí)是未知的，采用下面的預(yù)測(cè)-修正方法求得。

對(duì)于任意接觸點(diǎn)k，首先假設(shè)其在時(shí)間步l處為阻滯狀態(tài)，即2個(gè)接觸點(diǎn)之間沒(méi)有相對(duì)運(yùn)動(dòng)，因此kλx(l)=0，得到非線(xiàn)性力的預(yù)測(cè)值：

kλpre(l)=kλu(l)

(21)

再利用庫(kù)倫摩擦定律來(lái)判斷該點(diǎn)的真實(shí)接觸狀態(tài)，并用kλx(l)對(duì)非線(xiàn)性力進(jìn)行修正，共有如下3種情況：

1) 分離：|kλpre,N(l)|+|N0|≤0，其中N0代表施加在接觸點(diǎn)上的初始正壓力。此時(shí)接觸對(duì)之間不存在接觸力的作用，即kλ(l)=0，此時(shí)：

kλx(l)=kλu(l)

(22)

2) 阻滯：|kλpre,N(l)|+|N0|>0，同時(shí)切向力|kλpre,T(l)|<μf|kλpre,N(l)+N0|。其中μf代表接觸面的摩擦系數(shù)。滿(mǎn)足此條件時(shí)，說(shuō)明當(dāng)前的接觸狀態(tài)與預(yù)測(cè)結(jié)果相吻合，無(wú)需對(duì)非線(xiàn)性力進(jìn)行修正，即

kλx(l)=0

(23)

3) 滑移：|kλpre,N(l)|+|N0|>0，同時(shí)切向力|kλpre,T(l)|≥μf|kλpre,N(l)+N0|。此時(shí)法向接觸點(diǎn)之間相對(duì)位移的值為0；切向存在相對(duì)運(yùn)動(dòng)，且根據(jù)庫(kù)倫摩擦定律，切向摩擦力的大小可表示為μf|kλpre,N(l)+N0|，方向與相對(duì)速度的方向一致，因此可以得到滑移狀態(tài)下的修正項(xiàng)：

(24)

需要注意的是，基于牛頓-拉夫遜法的迭代方法為局部收斂法，即當(dāng)給的初值與方程的解較接近時(shí)，該方法具有較強(qiáng)的數(shù)值穩(wěn)定性，反之方程容易不收斂。因此在求解頻域內(nèi)的強(qiáng)迫響應(yīng)時(shí)，可以采用自適應(yīng)步長(zhǎng)的自然延拓法從小到大進(jìn)行掃頻計(jì)算，并以上一個(gè)頻率點(diǎn)的迭代收斂解作為計(jì)算下一個(gè)頻率點(diǎn)強(qiáng)迫響應(yīng)的迭代初值。對(duì)于第一個(gè)頻率點(diǎn)，一般將遠(yuǎn)離共振峰處系統(tǒng)不受非線(xiàn)性力作用下相對(duì)位移的頻域值作為迭代初值。

2.3 用解析雅克比矩陣加速計(jì)算

在利用牛頓-拉夫遜法以及其改進(jìn)方法進(jìn)行迭代時(shí)，其中必要的一步則是用非線(xiàn)性代數(shù)方程組式(15)的雅克比矩陣來(lái)判斷下一次迭代的值，在MATLAB中用fsolve函數(shù)解非線(xiàn)性方程組時(shí)，默認(rèn)情況下該矩陣一般通過(guò)有限差分法給出。然而當(dāng)方程中的非線(xiàn)性自由度數(shù)較大時(shí)，利用數(shù)值差分計(jì)算雅克比矩陣的時(shí)間成本較大，且精度也較低。如果能給出解析形式的雅克比矩陣，就能夠有效地解決這一問(wèn)題[27]。因此，給出適用于速度型動(dòng)態(tài)拉格朗日法的解析雅克比矩陣。

對(duì)于收斂殘差形式的運(yùn)動(dòng)方程(16)，其雅克比矩陣可寫(xiě)為

(25)

式中：Dr直接通過(guò)式(15)計(jì)算得到，因此計(jì)算J的關(guān)鍵就在于求解拉格朗日乘子對(duì)于相對(duì)位移矢量的偏導(dǎo)。利用求導(dǎo)的鏈?zhǔn)椒▌t，式(25)的第2項(xiàng)寫(xiě)為

(26)

其中：

(27)

圖3 基于時(shí)頻轉(zhuǎn)換的動(dòng)態(tài)拉格朗日法流程
Fig.3 Flowchart of DLFT based on time-frequency conversion

(28)

(29)

(30)

?kλx,N(l)/?hλu,T(l)的值為0。另外2種情況則為切向的修正項(xiàng)關(guān)于切向的預(yù)測(cè)項(xiàng)的偏導(dǎo)?kλx,T(l)/?hλu,T(l)：

(31)

以及切向的修正項(xiàng)關(guān)于法向的預(yù)測(cè)項(xiàng)的偏導(dǎo)?kλx,T(l)/?hλu,N(l)：

(32)

(33)

(34)

通過(guò)式(25)～式(34)即可獲得雅克比矩陣的解析表達(dá)式，從而提高利用DLFT法預(yù)測(cè)系統(tǒng)強(qiáng)迫響應(yīng)的計(jì)算效率與數(shù)值穩(wěn)定性。

2.4 計(jì)算流程簡(jiǎn)述

圖4給出了針對(duì)具有摩擦片的有限元模型強(qiáng)迫響應(yīng)預(yù)測(cè)算法的計(jì)算流程，概述如下：① 利用標(biāo)準(zhǔn)有限元軟件(本文使用ANSYS)對(duì)帶有波紋形干摩擦阻尼片的結(jié)構(gòu)的線(xiàn)性部分進(jìn)行建模；② 提取系統(tǒng)的剛度矩陣K、質(zhì)量矩陣M與阻尼矩陣C；③ 對(duì)得到的K、M、C矩陣進(jìn)行固定界面模態(tài)減縮，僅保留少量的主節(jié)點(diǎn)自由度。值得注意的是，該減縮僅在生成待求解代數(shù)方程組的過(guò)程中進(jìn)行一次；④ 將方程的自由度凝聚到非線(xiàn)

圖4 基于減縮模型的DLFT計(jì)算流程
Fig.4 Flowchart of DLFT based on reduced order model

性自由度上，進(jìn)一步減少計(jì)算規(guī)模；⑤ 利用圖3所示DFLT法求解非線(xiàn)性力的大小，在該步驟中采用時(shí)頻轉(zhuǎn)換技術(shù)，即先將頻域獲得的位移轉(zhuǎn)化到時(shí)域，并在時(shí)域內(nèi)對(duì)非線(xiàn)性力進(jìn)行修正，再將修正的非線(xiàn)性力轉(zhuǎn)回頻域；⑥ 將非線(xiàn)性力代回運(yùn)動(dòng)方程式(10)，得到觀察點(diǎn)的強(qiáng)迫響應(yīng)計(jì)算結(jié)果。

2.5 程序正確性測(cè)試

對(duì)于2.4節(jié)計(jì)算流程，可以通過(guò)如圖5所示的具有2個(gè)摩擦面的集中參數(shù)模型，將計(jì)算結(jié)果與基于切/法向剛度的時(shí)頻轉(zhuǎn)換技術(shù)(Elastic Frequency-Time technique, EFT)的結(jié)果進(jìn)行比較，從而說(shuō)明程序的正確性。

圖5中的模型分為上下2個(gè)結(jié)構(gòu)，不同結(jié)構(gòu)的質(zhì)量塊之間存在接觸。在驗(yàn)證模型中，每一個(gè)接觸點(diǎn)包含1個(gè)切向與1個(gè)法向。在上結(jié)構(gòu)中，每一個(gè)質(zhì)量塊受到大小為F的激振力激勵(lì)。模型中各參數(shù)的取值如表1所示。

圖5 帶有接觸面的集中參數(shù)模型
Fig.5 Lumped parameter model with contact interfaces

表1 帶有接觸的集中參數(shù)模型參數(shù)

利用切向/法剛度法驗(yàn)證速度型動(dòng)態(tài)拉格朗日算法的正確性時(shí)，用如圖6所示的考慮接觸剛度的遲滯庫(kù)倫摩擦模型來(lái)描述接觸面的力學(xué)特征。該模型認(rèn)為在阻滯狀態(tài)下，由于力的作用，接觸點(diǎn)之間仍存在一定的彈性變形，并通過(guò)引入切/法向彈簧來(lái)描述這一特征。在驗(yàn)證模型中，切向彈簧與法向彈簧的剛度分別為kt，kn=1×107N/m。

圖6 二維遲滯彈簧接觸模型示意圖
Fig.6 Illustration of 2D contact model with hysteresis springs

分別利用速度型動(dòng)態(tài)拉格朗日法與切/法向剛度法計(jì)算該模型在不同正壓力下的強(qiáng)迫響應(yīng)，其中觀測(cè)點(diǎn)為質(zhì)量塊1沿著x方向的位移，同時(shí)上述2種方法都是基于解析雅克比矩陣的。兩者的結(jié)果對(duì)比如圖7所示。

圖7中標(biāo)注的離散點(diǎn)為利用基于切/法向遲滯彈簧模型的NewMark法直接在時(shí)域內(nèi)求得的特定頻率下的響應(yīng)幅值。從圖7可以看出，利用切向/法向剛度法與速度型動(dòng)態(tài)拉格朗日法的計(jì)算結(jié)果十分相近，其中存在的微小的差異是由于DLFT在時(shí)域上是利用不考慮遲滯彈簧的庫(kù)倫模型去滿(mǎn)足接觸約束的。該計(jì)算結(jié)果與文獻(xiàn)[33]中描述的現(xiàn)象相吻合。在計(jì)算時(shí)間上，以圖7為例，利用速度型動(dòng)態(tài)拉格朗日法得到一條穩(wěn)態(tài)響應(yīng)曲線(xiàn)平均所用的CPU時(shí)間為145.3 s，而利用切/法向剛度法平均所用的CPU時(shí)間為638.7 s，可見(jiàn)在計(jì)算結(jié)果近似的情況下，前者的計(jì)算效率要明顯高于后者。

圖7 不同算法下強(qiáng)迫響應(yīng)曲線(xiàn)的比較
Fig.7 Comparison of forced response curves using different algorithms

3 波紋摩擦片的減振效果分析

本節(jié)將用上述數(shù)值工具研究波紋形摩擦阻尼片對(duì)于薄壁結(jié)構(gòu)的減振效果。注意到，干摩擦阻尼器的減振效果與系統(tǒng)在全阻滯與全滑移狀態(tài)下的共振頻率差相關(guān)。一般來(lái)說(shuō)較大的頻率差對(duì)應(yīng)著較好的減振性能[37]。因此阻尼片應(yīng)該布置在能夠明顯調(diào)整主體結(jié)構(gòu)在目標(biāo)頻帶內(nèi)的固有頻率的位置上。

平板結(jié)構(gòu)為一種具有代表性的薄壁結(jié)構(gòu)，本節(jié)通過(guò)在一塊405 mm×405 mm，厚度為2 mm的板上布置波紋形干摩擦阻尼片，利用仿真模擬的手段研究阻尼片減振性能。板由鋼構(gòu)成，其材料參數(shù)為：彈性模量E=210 GPa，泊松比μ=0.3，密度ρ=7 800 kg/m3。邊界條件為板的四邊全部固支。

將第5階模態(tài)作為目標(biāo)模態(tài)，研究阻尼片對(duì)于由此模態(tài)所主導(dǎo)的振動(dòng)的減振效果。其模態(tài)頻率為400.47 Hz，模態(tài)振型如圖8所示。將4片阻尼片布置到如圖9所示的位置，每片阻尼片都對(duì)應(yīng)平板四條邊的中點(diǎn)，且其中一個(gè)接觸面位于固定端3 mm處。阻尼片的材料參數(shù)為：彈性模量E=70 GPa，泊松比μ=0.3，密度ρ=2 700 kg/m3。

圖8 四邊固支板的第5階模態(tài)振型
Fig.8 Modal shape of the 5th mode of four side fixed plate

阻尼片的幾何尺寸如表2所示。為激起該階模態(tài)，在圖9中的施力點(diǎn)位置施加垂直于板方向的激振力，大小為F=15 N，同時(shí)相鄰2點(diǎn)之間的激振力存在π/2的相位差。

圖9 帶有4個(gè)阻尼片的板結(jié)構(gòu)
Fig.9 Plate with four frictional patches attached

在分析時(shí)，首先建立圖9中含有阻尼片的板結(jié)構(gòu)的有限元模型(利用ANSYS)；再用固定界面模態(tài)減縮來(lái)降低模型的自由度，保留其中施力點(diǎn)與觀測(cè)點(diǎn)的線(xiàn)性自由度，以及接觸面上的所有非線(xiàn)性自由度。為保證計(jì)算精度，截取前50階模態(tài)，由此得到的降階模型前10階的模態(tài)頻率與原模型相比，誤差小于0.5%，且模態(tài)振型相似性約等于1，故可以用此降階模型的計(jì)算結(jié)果來(lái)替代原模型。

利用DLFT方法對(duì)上述模型進(jìn)行計(jì)算，模型本身的阻尼比取ξ=0.002，文獻(xiàn)[38]中表明摩擦系數(shù)的大小對(duì)干摩擦阻尼器的最優(yōu)減振效果并不敏感，這里該參數(shù)的取值為μf=0.6。為保證計(jì)算精度，本文中保留的諧波數(shù)為Nh=5，計(jì)算得到如圖10所示的不同正壓力下觀測(cè)點(diǎn)的穩(wěn)態(tài)強(qiáng)迫響應(yīng)曲線(xiàn)。隨著正壓力的增大，共振頻率向右偏移，完全自由狀態(tài)下的共振頻率為400.6 Hz，與純平板結(jié)構(gòu)下的固有頻率變化不大，說(shuō)明當(dāng)阻尼片布置在該位置時(shí)，其附加質(zhì)量對(duì)整體結(jié)構(gòu)的動(dòng)力學(xué)特性幾乎沒(méi)有影響。當(dāng)接觸面處于完全阻滯狀態(tài)時(shí)，共振頻率為408.4 Hz，頻率偏移量為7.8 Hz。共振幅值隨著正壓力的增加先減小，后增大。當(dāng)初始正壓力為80 N時(shí)，阻尼片的減振性能最佳，共振幅值相對(duì)于完全分離狀態(tài)下降了75.4%。同時(shí)，阻尼片的質(zhì)量(考慮了安裝所用的螺栓和螺母)占主體結(jié)構(gòu)的0.6%。由此可見(jiàn)，該波紋形干摩擦阻尼器對(duì)于平板結(jié)構(gòu)具有較為理想的減振性能。

還可以對(duì)接觸區(qū)域的局部動(dòng)力學(xué)特性進(jìn)行分析。以圖9中箭頭所指的接觸點(diǎn)為例，在初始正壓力為80 N下達(dá)到共振峰時(shí)，該接觸點(diǎn)在時(shí)域內(nèi)的運(yùn)動(dòng)狀態(tài)如圖11所示。其中實(shí)線(xiàn)代表接觸點(diǎn)處所受的切向摩擦力，虛線(xiàn)代表接觸對(duì)之間的切向相對(duì)位移。當(dāng)接觸狀態(tài)為阻滯時(shí)(陰影部分)，接觸對(duì)之間的相對(duì)位移保持不變；當(dāng)接觸狀態(tài)為滑移時(shí)，由于切法向運(yùn)動(dòng)耦合的影響，導(dǎo)致切向的摩擦力并不是定值。

關(guān)于計(jì)算效率，以最優(yōu)正壓力下計(jì)算強(qiáng)迫響應(yīng)為例，所用的時(shí)間與減縮過(guò)程中系統(tǒng)的自由度的變化如表3所示。可見(jiàn)構(gòu)造減縮模型所耗費(fèi)的時(shí)間只占總分析時(shí)間的很小部分，卻能有效地降低模型的自由度(降低了99.7%)，因此推薦在用DLFT分析較大規(guī)模的有限元模型時(shí)使用減縮技術(shù)。如果不采用解析雅克比矩陣而是用默認(rèn)的數(shù)值差分近似，平均一個(gè)頻率點(diǎn)所需的CPU時(shí)間為333.7 s；采用減縮方法后，計(jì)算706個(gè)頻率點(diǎn)總共花了5 360.1 s，平均每1個(gè)頻率點(diǎn)的CPU計(jì)算時(shí)間為7.6 s?？梢?jiàn)如式(25)～式(34)所示的解析雅克比矩陣對(duì)速度型的動(dòng)態(tài)拉格朗日法的計(jì)算效率至關(guān)重要，在計(jì)算自由度相對(duì)較多的結(jié)構(gòu)時(shí)幾乎是不可或缺的。

圖10 不同正壓力下的穩(wěn)態(tài)強(qiáng)迫響應(yīng)結(jié)果
Fig.10 Steady state response of system under different normal preloads

圖11 接觸點(diǎn)切向力和相對(duì)位移的時(shí)間歷程
Fig.11 Time-histories of tangential force and relative displacement for contact node

表3 計(jì)算時(shí)間與減縮過(guò)程中自由度數(shù)的變化

4 結(jié) 論

本文提出了一種適用于薄壁結(jié)構(gòu)的波紋形干摩擦阻尼器，其具有結(jié)構(gòu)簡(jiǎn)單、附加質(zhì)量小、易于安裝與調(diào)節(jié)初始正壓力的優(yōu)點(diǎn)。利用仿真模擬的手段對(duì)阻尼片的減振性能進(jìn)行了分析，在高階諧波平衡法的基礎(chǔ)上引入了速度型動(dòng)態(tài)拉格朗日法，并通過(guò)提供解析形式的雅克比矩陣來(lái)提高迭代效率與計(jì)算穩(wěn)定性，從而實(shí)現(xiàn)了對(duì)帶有接觸系統(tǒng)的頻域強(qiáng)迫響應(yīng)的快速預(yù)測(cè)，之后的仿真過(guò)程表明該算法具有較高的計(jì)算效率。

利用強(qiáng)迫響應(yīng)預(yù)測(cè)算法對(duì)帶有阻尼片的薄壁結(jié)構(gòu)進(jìn)行分析，可以看出，對(duì)于平板而言，通過(guò)合理布置阻尼片的位置，增大結(jié)構(gòu)自由/阻滯時(shí)的頻率差，從而獲得較好的減振效果。以文中的平板結(jié)構(gòu)為例，僅使用質(zhì)量占整體結(jié)構(gòu)0.6%的阻尼片就能使得共振峰值下降75.4%。

現(xiàn)有工作的一個(gè)展望是，針對(duì)干摩擦阻尼器在各頻段下最優(yōu)正壓力不同的特點(diǎn)，還可以引入主動(dòng)控制環(huán)節(jié)通過(guò)設(shè)計(jì)可控正壓力提高干摩擦阻尼器在寬頻內(nèi)的減振性能[30，39]。這些后續(xù)工作仍然需要一種可以考慮干摩擦界面切-法向耦合的高效數(shù)值工具，本文提出的基于減縮模型和解析雅克比矩陣的速度型動(dòng)態(tài)拉格朗日乘子方法仍適用。

附錄A

布爾矩陣是元素只取0或1的矩陣，因此又被稱(chēng)為0-1矩陣。它的作用是提取含有m個(gè)元素的向量X中的某些特定項(xiàng)，并將提取的元素重新組成維度為n的向量Y。該過(guò)程所對(duì)應(yīng)的布爾矩陣B的維度為n×m。

記R為X和Y上的對(duì)應(yīng)關(guān)系，則布爾矩陣中B的各元素滿(mǎn)足：

圖A1X與Y之間的對(duì)應(yīng)關(guān)系
Fig.A1 Correspondences betweenXandY

(A1)

其轉(zhuǎn)置矩陣BT則能將Y中的元素放回X中的對(duì)應(yīng)位置。