基于掃頻響應(yīng)反推黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型

孫 偉， 王 茁， 朱明偉

(東北大學(xué)機械工程與自動化學(xué)院 沈陽，110819)

引 言

利用黏彈性材料耗能實現(xiàn)對結(jié)構(gòu)件的振動抑制已成為一種重要的被動阻尼減振方法。創(chuàng)建可以有效預(yù)估黏彈性復(fù)合結(jié)構(gòu)動力學(xué)特性的模型是該阻尼減振技術(shù)的基礎(chǔ)，要獲得分析模型，必須首先知道黏彈性材料的力學(xué)特性參數(shù)。黏彈性阻尼材料的力學(xué)特性參數(shù)包括儲能模量和損耗因子(或耗能模量)，具有頻率依賴性的特點，這是黏彈性阻尼材料最獨特的力學(xué)特點。黏彈性材料的頻率依賴性使有效辨識材料的力學(xué)特性參數(shù)成為一項具有挑戰(zhàn)性的研究任務(wù)。

目前，研究人員已提出大量的模型來表征黏彈性材料的頻率依賴性，典型的模型包括：復(fù)模量模型[1]、Golla-Hughes-McTavish模型[2]、滯彈性位移場模型[3]和分?jǐn)?shù)階導(dǎo)數(shù)模型[4]等。這使辨識黏彈性材料參數(shù)變?yōu)楂@得上述表征模型的各個特性參數(shù)。在上述模型中，分?jǐn)?shù)階導(dǎo)數(shù)模型可以有效表征真實黏彈性材料在時域和頻域的力學(xué)行為，因而得到了廣泛應(yīng)用[5-7]。

將數(shù)值分析和振動實驗相結(jié)合的反推法是辨識黏彈性材料力學(xué)特性參數(shù)的一種常用方法，很多學(xué)者針對黏彈性復(fù)合梁結(jié)構(gòu)進(jìn)行了反推辨識。Barkanov等[8]基于響應(yīng)面優(yōu)化技術(shù)，通過實測復(fù)合梁的振動衰減響應(yīng)來辨識用復(fù)模量模型表征的黏彈性材料參數(shù)。Shi等[9]通過實測復(fù)合梁的共振頻率和模態(tài)損耗因子，利用數(shù)值優(yōu)化法反推黏彈性材料的彈性模量和損耗因子。上述研究最大的缺陷是僅能針對共振頻率處的數(shù)據(jù)點進(jìn)行參數(shù)辨識，而由于數(shù)據(jù)點數(shù)量上的欠缺可能導(dǎo)致辨識的黏彈性參數(shù)與真實值存在偏差。文獻(xiàn)[10]針對一個自由阻尼層懸臂梁，用力錘激勵測得頻響函數(shù)，采用兩步法反推出用分?jǐn)?shù)階導(dǎo)數(shù)模型表征的黏彈性材料參數(shù)。該方法已經(jīng)不局限于僅針對共振頻率點進(jìn)行參數(shù)辨識，是黏彈性材料力學(xué)特性參數(shù)辨識方法上的一個巨大進(jìn)步。但是，錘擊測試涉及激勵位置難以準(zhǔn)確確定的問題，不精確的激勵點可能影響?zhàn)椥圆牧蠀?shù)的辨識結(jié)果。另外，上述研究在創(chuàng)建的理論模型中只考慮了單一的黏彈性材料阻尼，沒有考慮梁夾持端以及空氣中存在的阻尼(這里定義為剩余等效黏性阻尼)。這樣的模型將導(dǎo)致測量的阻尼值總是大于計算的阻尼值，將影響?zhàn)椥圆牧蠀?shù)的辨識結(jié)果。

筆者提出一種基于掃頻響應(yīng)來反推黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)的方法。該方法具有以下特點：a. 以具有自由阻尼層的懸臂板形試件為對象開展黏彈性材料參數(shù)辨識研究；b. 理論模型中同時考慮了材料阻尼和剩余等效黏性阻尼；c. 利用共振點及非共振點處的數(shù)據(jù)實施分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)辨識；d. 利用基礎(chǔ)激勵下掃頻得到的頻域響應(yīng)進(jìn)行參數(shù)辨識。研究中建立了基礎(chǔ)激勵作用下黏彈性復(fù)合板有限元方程，給出了求解振動響應(yīng)的方法。提出面向振動響應(yīng)、基于靈敏度法的匹配計算來實現(xiàn)黏彈性參數(shù)的反推辨識。最后，以貼敷ZN-1型黏彈性材料的懸臂鈦板為對象，實現(xiàn)了該黏彈性材料用分?jǐn)?shù)階導(dǎo)數(shù)模型描述的力學(xué)特性參數(shù)辨識。

1 反推法辨識原理

如圖1所示，基于反推法的黏彈性材料參數(shù)辨識主要涉及3項技術(shù)，分別為振動測試、理論建模及匹配計算。待辨識的黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型可表示為

(1)

圖1 基于反推法辨識黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型的原理Fig.1 The procedure of identifying the fractional derivative model of viscoelastic materials

對于不同的聚合物組合以及承受的力學(xué)狀態(tài)的差異，通常黏彈性材料表現(xiàn)出不同的力學(xué)特點[10]。這里選擇貼敷黏彈性材料的板型試件來辨識黏彈性材料參數(shù)，是因為基于梁型試件辨識的結(jié)果不一定完全滿足板型試件動力學(xué)分析的需要。

辨識原理中的振動測試包含兩部分：a.獲得阻尼處理前懸臂板的各階模態(tài)阻尼比；b.獲得基礎(chǔ)激勵作用下黏彈性阻尼板的頻域掃頻響應(yīng)。其中，模態(tài)阻尼比將作為剩余等效黏性阻尼引入分析模型，而測得的頻域掃頻響應(yīng)主要用于反推辨識黏彈性材料參數(shù)。理論分析主要用于求解黏彈性阻尼板的頻域振動響應(yīng)，筆者創(chuàng)建的模型中同時考慮了黏彈性材料阻尼以及剩余等效黏彈性阻尼。

匹配計算是為了使理論模型獲得的振動響應(yīng)值與實驗值的偏差最小，進(jìn)而反推出所要獲得的分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)。這里提出了以振動響應(yīng)為目標(biāo)函數(shù)、基于靈敏度的匹配計算方法。需要說明的是，匹配計算可同時基于共振及非共振點進(jìn)行材料參數(shù)辨識，但是考慮到共振點信噪比較高，因而應(yīng)優(yōu)先選擇共振點。此外，選擇的共振及非共振頻率點總數(shù)應(yīng)大于分?jǐn)?shù)階導(dǎo)數(shù)模型待辨識參數(shù)的個數(shù)，對于本實例選擇的頻率點數(shù)應(yīng)大于4。

2 黏彈性復(fù)合板振動響應(yīng)分析

2.1 運動方程的建立及響應(yīng)求解

如圖2所示，當(dāng)薄板產(chǎn)生彎曲振動時，自由阻尼層會隨基體一起運動，在阻尼層內(nèi)部將產(chǎn)生拉壓變形而產(chǎn)生耗能，從而起到減振的作用。筆者認(rèn)為對于這樣的結(jié)構(gòu)，層間剪切效應(yīng)可以被忽略，而將黏彈性阻尼層與基體視為一體。事實上，文獻(xiàn)[10-11]也是按照這種忽略層間效應(yīng)來對黏彈性自由阻尼層復(fù)合結(jié)構(gòu)進(jìn)行建模的。

設(shè)坐標(biāo)軸恰好位于復(fù)合板中性面上，Hv為黏彈性阻尼層厚度，He為金屬基體厚度，Hev為基體及黏彈性層中面之間的距離，d為金屬基體中面到復(fù)合結(jié)構(gòu)中性面的距離。

圖2 基礎(chǔ)激勵下的黏彈性復(fù)合懸臂板Fig.2 The cantilever plate attached with viscoelastic FLD under base excitation

對式(1)的分?jǐn)?shù)階導(dǎo)數(shù)模型進(jìn)行實部及虛部分解，可得到復(fù)模量模型，其表達(dá)式為

(2)

用此復(fù)模量模型進(jìn)行黏彈性阻尼板的振動特性分析。由于在辨識算法中僅涉及若干頻率點對應(yīng)的振動響應(yīng)值，頻率信息已包含在內(nèi)，因而在響應(yīng)的求解過程中不需要考慮黏彈性材料的頻率依賴性。

對于處于基礎(chǔ)激勵作用下的該黏彈性復(fù)合板，其頻域運動方程為

(-ω2M+iωC+K*)X=F

(3)

其中：M為復(fù)合板的質(zhì)量矩陣；C為由夾持邊界及空氣產(chǎn)生的剩余等效黏性阻尼矩陣；K*為基體及黏彈性阻尼層共同產(chǎn)生的復(fù)剛度矩陣(包含了材料阻尼)；X和F分別為相對響應(yīng)向量和激振力向量。

X和F表達(dá)式為

X=U0-GU

(4a)

(4b)

這里采用模態(tài)疊加法求解黏彈性復(fù)合板在基礎(chǔ)激勵作用下的振動響應(yīng)。求解實模態(tài)的特征方程可表示為

(5)

求得的各階模態(tài)振型向量可組成模態(tài)振型矩陣φ，用模態(tài)振型矩陣φ對式(3)進(jìn)行解耦，得到一組獨立的以模態(tài)坐標(biāo)xNr(r=1,2,…,n)表達(dá)的運動方程，即

(6)

式(6)括號內(nèi)第1項和第3項系數(shù)分別為

其中：ηr為第r階模態(tài)損耗因子，來自于黏彈性材料阻尼。

根據(jù)振型阻尼假設(shè)，式(6)中括號內(nèi)第2項可以寫為

(8)

其中：ξr為第r階剩余等效模態(tài)阻尼比，可按照結(jié)構(gòu)阻尼處理前的各階模態(tài)阻尼比引入。

將式(7)，(8)代入式(6)，可以將方程轉(zhuǎn)換到模態(tài)坐標(biāo)下，整理后得到

(9)

由式(9)可知，在響應(yīng)的求解中同時考慮了兩種阻尼：a. 由黏彈性材料產(chǎn)生的模態(tài)損耗因子ηr；b. 由剩余等效黏性阻尼產(chǎn)生的模態(tài)阻尼比ξr。

(10)

結(jié)構(gòu)在基礎(chǔ)激勵作用下振動響應(yīng)為

(11)

其中：|·|為求模運算。

實際計算時一般不用取所有n階模態(tài)，考慮的階次m只需大于分析頻率范圍內(nèi)結(jié)構(gòu)的最高階次即可。

2.2 剛度及質(zhì)量矩陣的生成

對于圖2所示的黏彈性自由阻尼復(fù)合板，選用四邊形板單元(該單元共有4個節(jié)點，每個節(jié)點有3個自由度)進(jìn)行模擬。按照以下流程求解剛度矩陣和質(zhì)量矩陣。

1) 將薄板復(fù)合結(jié)構(gòu)劃分為若干個單元。

2) 求解單元的剛度及質(zhì)量矩陣，求解式分別為

De和Dv分別為基體和黏彈性阻尼層的彈性矩陣，表達(dá)式為

(13a)

(13b)

其中：μe和μv分別為彈性基體和黏彈性層的泊松比。

復(fù)合板的中性面到基體中心的距離d表示為

d=EvHvHev/(EeHe+EvHv)

(14)

3) 在考慮約束條件的基礎(chǔ)上組合單元剛度矩陣及單元質(zhì)量矩陣，最終可形成黏彈性復(fù)合板的總剛度矩陣K*和總質(zhì)量矩陣M。

3 基于靈敏度的匹配計算

參照式(1)，將待辨識的分?jǐn)?shù)階導(dǎo)數(shù)模型的4個參數(shù)E0，E∞，τ和α作為設(shè)計變量，采用差分法獲得靈敏度。以E0為例，第i個頻率點振動響應(yīng)對設(shè)計變量E0的靈敏度SE0,i可表示為

(15)

其中：ΔE0為設(shè)計變量的變化量(即步長)，這里設(shè)為10-4E0；ΔxE0,i為設(shè)計變量變化ΔE0時振動響應(yīng)的變化量(殘差)。

同樣，分別求出振動響應(yīng)對E∞，τ和α的靈敏度SE∞,i，Sτ,i和Sα,i。假設(shè)在整個辨識中共選擇P個頻率點，則求得的針對設(shè)計變量的靈敏度可組成靈敏度矩陣S，即

(16)

該匹配計算的目標(biāo)函數(shù)為

(17)

其中：Wx，Wb分別為振動響應(yīng)和設(shè)計變量的權(quán)重矩陣；Δb為設(shè)計變量變化量所組成的向量；Δx為對應(yīng)各頻率點振動響應(yīng)的殘差組成的向量。

Δb和Δx的表達(dá)式為

(18a)

(18b)

需要說明的是，式(18)中的殘差對應(yīng)所有設(shè)計變量變化時響應(yīng)的變化量。

得到如下迭代公式

(19)

經(jīng)過若干次迭代，滿足收斂條件后迭代終止，獲得對應(yīng)分?jǐn)?shù)階導(dǎo)數(shù)模型待辨識的4個參數(shù)。需要說明的是，考慮到共振點對應(yīng)的響應(yīng)值信噪比高，設(shè)置權(quán)重時應(yīng)將共振點的權(quán)重設(shè)為大值，將非共振點權(quán)重設(shè)為小值。

4 實例研究

以貼敷ZN-1型黏彈性材料的懸臂鈦板為例，基于掃頻獲得頻域響應(yīng)反推分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)。

4.1 試件及振動特性測試

阻尼處理前后的鈦板如圖3所示。該薄板的長、寬分別為154.1 mm，109.6 mm，鈦板基體的厚度為1.32 m，薄板夾持區(qū)長度為20 mm，擰緊力矩為32 N·m。在薄板的一側(cè)貼敷ZN-1型黏彈性材料，厚度為0.5 mm。鈦板的材料參數(shù)是已知的：儲能模量為110.32 GPa，密度為4 420 kg/m3，損耗因子為0.000 7。ZN-1型黏彈性材料的密度可實測獲得，具體為789.5 kg/m3，黏彈性材料及鈦板的泊松比分別取0.49和0.3，ZN-1型黏彈性材料的分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)需要辨識獲得。

實驗現(xiàn)場如圖4所示，涉及的儀器設(shè)備如表1所示。整個測試過程中均用激光測振儀拾振，參照圖4中的坐標(biāo)，拾振點位置為x=11 mm，y=67 mm。在基于反推法的辨識中，實驗系統(tǒng)與理論模型相對應(yīng)是至關(guān)重要的。剔除夾持區(qū)后，實驗結(jié)構(gòu)為受振動臺基礎(chǔ)激勵的、單面全部貼敷ZN-1型阻尼材料的和由夾具固定的(實踐表明，32 N·m的螺栓擰緊力矩可以保證板的底部為完全固定約束)懸臂結(jié)構(gòu)系統(tǒng)，這就保證了實驗系統(tǒng)與理論模型之間激勵條件、結(jié)構(gòu)形式和邊界條件均具有一致性。

圖3 鈦板實驗件Fig.3 Specimens used in the experiment

圖4 測試現(xiàn)場圖Fig.4 Actual photograph of vibration test

表1 主要測試儀器Tab.1 The instruments used in this test

對阻尼處理前的薄板進(jìn)行錘擊測試，采用半功率帶寬法由獲得的頻響函數(shù)辨識出各階模態(tài)阻尼比，測試結(jié)果如表2所示。利用振動臺，設(shè)定包含共振頻率的掃頻區(qū)間，對阻尼處理后的薄板進(jìn)行分段掃頻測試。選擇第1，3，5階共振區(qū)的頻域響應(yīng)用于黏彈性材料參數(shù)辨識，對應(yīng)的激勵幅度分別為0.25，1和2 g，相應(yīng)的測試結(jié)果如圖5所示。

表2阻尼處理前薄板的模態(tài)阻尼比

Tab.2Modaldampingratiosofthetitaniumplatewithoutdampingtreatment%

4.2 黏彈性材料的力學(xué)特性參數(shù)辨識

參照圖5所示的掃頻響應(yīng)，選擇3個共振頻率點以及每個共振頻率左右各一個非共振頻率點(即3個共振點，6個非共振點)，利用本研究方法進(jìn)行參數(shù)辨識。

圖5 實驗測試及理論計算得到的頻域響應(yīng)Fig.5 Frequency domain responses obtained by experiment and analysis

整個模型共劃分了100個單元，121個節(jié)點，利用Matlab編程完成整個計算。參考文獻(xiàn)[13]選擇ZN-1型黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型的初值，具體為：E0=1.8 MPa，E∞=12.5 MPa，τ=2.42×10-4，α=0.99。計算各頻率點處振動響應(yīng)對分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)的靈敏度以及振動響應(yīng)的殘差向量。接著，將分?jǐn)?shù)階導(dǎo)數(shù)模型各參數(shù)的權(quán)重均取0.25，在振動響應(yīng)的權(quán)重矩陣中，共振點處的權(quán)重取0.2，非共振點處的權(quán)重取0.06。最后，將確定的靈敏度矩陣、振動響應(yīng)殘差向量以及權(quán)重矩陣代入式(17)進(jìn)行迭代計算。設(shè)收斂條件為：計算得到的振動響應(yīng)與相應(yīng)的實驗值偏差小于1%。經(jīng)過若干次迭代，可反推獲得ZN-1型黏彈性材料的分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)，具體結(jié)果如表3所示。例如，以激振頻率為橫軸， ZN-1型黏彈性材料的儲能模量或損耗因子為縱軸，則可對該材料進(jìn)行頻率依賴性頻域表征，如圖6所示。

表3ZN-1型黏彈性材料的分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)

Tab.3FractionalderivativemodelparametersofZN-1viscoelasticmaterial

參數(shù)E0E∞τα辨識值1．5MPa11．3MPa1．02×10-40．97

圖6 ZN-1型黏彈性材料的頻率依賴性表征Fig.6 Characterization of ZN-1 viscoelastic material with frequency dependent characteristics

為了說明辨識結(jié)果的合理性，將獲得的黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型寫成復(fù)模量的形式并代入到黏彈性阻尼板振動響應(yīng)分析模型中。黏彈性復(fù)合板的頻域振動響應(yīng)需按照各頻率點進(jìn)行求解(計算步長為1 Hz)。首先，針對每個頻率值確定黏彈性材料的儲能模量和損耗因子；其次，確定總剛度矩陣K*和總質(zhì)量矩陣M，并進(jìn)行模態(tài)分析；然后，取8階模態(tài)，將表2中的阻尼值作為剩余等效黏性阻尼輸入分析模型，按照模態(tài)疊加法求解每個頻率點對應(yīng)的響應(yīng)值；最后，依次類推求解各頻率點對應(yīng)的響應(yīng)值，繪制與實驗頻率范圍一致的頻域響應(yīng)曲線，如圖5所示。提取各指定頻率點對應(yīng)的理論計算與實驗測試獲得的振動響應(yīng)值并進(jìn)行比對，如表4所示?？梢钥闯?，理論計算獲得的黏彈性阻尼板頻域響應(yīng)值與實測值幾乎一致，證明了辨識結(jié)果的合理性。

表4實驗測試及理論計算獲得的振動響應(yīng)

Tab.4Vibrationresponsesobtainedbyexperimentandanalyticalcalculation

頻率點/Hz振動響應(yīng)/(m·s-1)實驗測試?yán)碚撚嬎闫?%56．540．03090．03041．6259．13?0．23510．23560．2162．030．02440．02471．23363．090．13080．13291．61367．11?0．48540．48400．28372．210．12350．12591．94642．100．01930．01993．11645．08?0．06440．06480．62649．360．01300．01353．85

*號為共振點

5 結(jié) 論

1) 在同時考慮黏彈性材料阻尼以及剩余等效黏性阻尼的基礎(chǔ)上，創(chuàng)建了黏彈性阻尼板振動響應(yīng)分析模型。以理論計算獲得的振動響應(yīng)與實測值偏差最小為目標(biāo)函數(shù)，研發(fā)出基于靈敏度的匹配算法，在同時保證計算精度和計算效率的前提下實現(xiàn)了黏彈性材料分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)的辨識。

2) 以ZN-1型黏彈性材料為例，用本研究方法辨識出該材料的分?jǐn)?shù)階導(dǎo)數(shù)模型參數(shù)。將獲得的辨識結(jié)果代入到黏彈性阻尼板振動響應(yīng)分析模型中，并比較計算得到的響應(yīng)值與實測值，發(fā)現(xiàn)兩者幾乎一致，從而證明了辨識結(jié)果的合理性。利用實例中9個頻率點(3個共振點，6個非共振點)辨識出的黏彈性材料參數(shù)對復(fù)合板結(jié)構(gòu)進(jìn)行動力學(xué)分析，完全可以再現(xiàn)實際的黏彈性阻尼板頻域振動響應(yīng)行為。

[1] Martinez-Agirre M, Elejabarrieta M J. Dynamic characterization of high damping viscoelastic materials from vibration test data[J]. Journal of Sound and Vibration, 2011, 330(16): 3930-3943.

[2] McTavish D J, Hughes P C. Modeling of linear viscoelastic space structures[J]. Journal of Vibration and Acoustics, 1993, 115(1): 103-110.

[3] Lesieutre G A, Bianchini E. Time domain modeling of linear viscoelasticity using anelastic displacement fields[J]. Journal of Vibration and Acoustics, 1995, 117(4): 424-430.

[4] Bagley R L, Torvik P J. Fractional calculus in the transient analysis of viscoelastically damped structures[J]. AIAA Journal, 1985, 23(6): 918-925.

[5] Escobedo-Torres J, Ricles J M. The fractional order elastic-viscoelastic equations of motion: formulation and solution methods[J]. Journal of Intelligent Material Systems and Structures, 1998, 9(7): 489-502.

[6] Schmidt A, Gaul L. Finite element formulation of viscoelastic constitutive equations using fractional time derivatives[J]. Nonlinear Dynamics, 2002, 29(1-4): 37-55.

[7] Cortés F, Elejabarrieta M J. Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models[J]. Computational Mechanics, 2007, 40(2): 313-324.

[8] Barkanov E, Skukis E, Petitjean B. Characterisation of viscoelastic layers in sandwich panels via an inverse technique[J]. Journal of Sound and Vibration, 2009, 327(3): 402-412.

[9] Shi Y, Sol H, Hua H. Material parameter identification of sandwich beams by an inverse method[J]. Journal of Sound and Vibration, 2006, 290(3): 1234-1255.

[10] Kim S Y, Lee D H. Identification of fractional-derivative-model parameters of viscoelastic materials from measured FRFs[J]. Journal of Sound and Vibration, 2009, 324(3): 570-586.

[11] Parthasarathy G, Reddy C V R, Ganesan N. Partial coverage of rectangular plates by unconstrained layer damping treatments[J]. Journal of Sound and Vibration, 1985, 102(2): 203-216.

[12] 劉建濤, 杜平安, 黃明鏡, 等. 阻尼連續(xù)體簡諧基礎(chǔ)振動的有限元實現(xiàn)方法[J]. 機械工程學(xué)報, 2010, 46(1): 109-114.

Liu Jiantao, Du Ping′an, Huang Mingjing, et al. Finite element approach for analyzing continua with damp under harmonic loads from the base[J]. Journal of Mechanical Engineering, 2010, 46(1): 109-114.(in Chinese)

[13] 李恩奇, 盛秀成, 王江楓. ZN-1型黏彈性阻尼材料模型參數(shù)修正研究[J]. 強度與環(huán)境, 2011, 38(5): 36-41.

Li Enqi, Sheng Xiucheng, Wang Jiangfeng. Parameter update of ZN-1 viscoelastic damping material model[J]. Structure & Environment Engineering, 2011, 38(5): 36-41. (in Chinese)