高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析的小波有限元方法實(shí)現(xiàn)

耿 佳,李 明,張興武,楊來浩,陳雪峰

(1.西安交通大學(xué) 機(jī)械工程學(xué)院,西安 710049；2.西安交通大學(xué) 機(jī)械制造系統(tǒng)工程國家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

數(shù)值分析方法因其具有可操作性強(qiáng)、成本低、可突破各種理論推導(dǎo)和實(shí)驗(yàn)研究條件限制等優(yōu)點(diǎn)，已被廣泛應(yīng)用于工程領(lǐng)域中進(jìn)行結(jié)構(gòu)振動(dòng)特征分析[1-4]。然而，近年來寬頻振動(dòng)分析問題已經(jīng)成為限制高端裝備發(fā)展的主要障礙[5-6]。著名學(xué)者M(jìn)ace和Desmet在聲振分析領(lǐng)域內(nèi)頂尖雜志Journal of Sound and Vibration中聯(lián)合聲明，在潛艇、火箭和飛機(jī)等高端裝備中包含有大量的薄殼、薄板和曲殼結(jié)構(gòu)，在該領(lǐng)域內(nèi)將此類結(jié)構(gòu)特征尺寸d遠(yuǎn)大于厚度t(即d/t>10)的結(jié)構(gòu)定義為高模態(tài)密度結(jié)構(gòu)，而在對上述結(jié)構(gòu)振動(dòng)特征進(jìn)行分析時(shí)普遍存在寬頻振動(dòng)分析問題[7-8]。

該問題存在的主要原因之一，是以傳統(tǒng)有限元方法為代表的確定性分析方法，在對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析時(shí)由于計(jì)算成本過高，耗散誤差明顯和參數(shù)不確定性等原因，而存在高頻域難以提供有效數(shù)值解的問題，如圖1左側(cè)所示，其中fTFEM為傳統(tǒng)有限元方法有效頻域上限。而以統(tǒng)計(jì)能量分析(statistic energy analysis，SEA)方法為代表的不確定性分析方法由于在帶寬內(nèi)的模態(tài)數(shù)大于5時(shí)才能提供有效數(shù)值解而難以在較低頻域內(nèi)進(jìn)行有效地振動(dòng)分析，如圖1右側(cè)所示，其中fSEA為不確定性分析方法有效頻域下限。并且，在fTFEM和fSEA之間存在一個(gè)頻率區(qū)域，正如圖中灰色區(qū)域所示，該頻域即為中頻域Ωmid，其可寫為

Ωmid=fSEA-fTFEM

(1)

顯然，確定性分析方法和不確定性分析方法均難以在該頻域內(nèi)對高模態(tài)密度結(jié)構(gòu)進(jìn)行有效地振動(dòng)分析，導(dǎo)致難以實(shí)現(xiàn)寬頻振動(dòng)分析[9]。

圖1 高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析Fig.1 Dynamic analysis in the wide-frequency domain

為了解決該問題，基于數(shù)值分析方法對高模態(tài)密度結(jié)構(gòu)實(shí)現(xiàn)有效的寬頻振動(dòng)分析。以魯汶大學(xué)國際聲振分析領(lǐng)軍人物Wim Desmet教授為代表的學(xué)者們，提出了更高求解效率的研究思路[10]。據(jù)此，為了提高有限元分析方法在更高頻域內(nèi)的求解精度，Harari等提出了h-FEA和p-FEA方法，用于拓寬有限元分析方法在進(jìn)行振動(dòng)分析時(shí)，能夠提供有效數(shù)值解的頻域范圍，其中h-FEA是通過細(xì)化網(wǎng)格來獲取更高的求解精度[11]，p-FEA是在特定網(wǎng)格劃分條件下使用更高階次多項(xiàng)式代替低階插值多項(xiàng)式的方式提高有限元求解精度[12]。上述兩種方法可在較低頻域內(nèi)緩解耗散誤差導(dǎo)致的數(shù)值解偏差較大的問題。Dai等提出了光滑有限元法，并進(jìn)行了相關(guān)應(yīng)用研究，弱化了耗散誤差帶來的影響[13]。隨后，研究者基于光滑有限元理論提出了節(jié)點(diǎn)光滑有限元法，邊界光滑有限元法等，并將其應(yīng)用于更復(fù)雜的振動(dòng)分析問題中，有效提升了有限元分析方法在進(jìn)行振動(dòng)分析時(shí)的頻域上限。Chazot等引入單元分解有限元方法，并結(jié)合平面波函數(shù)對短波振動(dòng)分析方程進(jìn)行了求解，在一定程度上提高了確定性分析方法的分析頻域上限[14]。

與此同時(shí)，為了緩解數(shù)值模型存在參數(shù)不確定性，導(dǎo)致有限元分析方法無法在高頻域內(nèi)提供有效數(shù)值解的問題，Chen等提出了將變量視為常數(shù)，引入隨機(jī)矩陣和向量，并通過一階泰勒展開進(jìn)行研究，基于攝動(dòng)理論將該泰勒展開式引入至有限元分析方法中進(jìn)行振動(dòng)分析研究[15]。Li等提出了將不確定性參數(shù)分析方法引入至光滑有限元的振動(dòng)分析方法，在一定程度上可在工程結(jié)構(gòu)存在參數(shù)不確定時(shí)提供有效的振動(dòng)分析數(shù)值解[16]。同時(shí)，為了更進(jìn)一步獲得具有隨機(jī)性的有限元建模方法，研究者們提出了隨機(jī)有限元方法，其在近些年受到了廣泛的關(guān)注，可以基于統(tǒng)計(jì)輸入?yún)?shù)來解決輸入?yún)?shù)不確定時(shí)遇到的數(shù)值解存在較大誤差的問題。

可以看出，經(jīng)過國內(nèi)外學(xué)者們長期以來的探索，無論是從降低耗散誤差，提高求解效率，還是不確定性參數(shù)建模方向，都獲得了豐富的研究成果，用以提升振動(dòng)分析過程中的頻域上限。然而，由于傳統(tǒng)有限元分析方法的理論限制，導(dǎo)致其無法大幅度降低耗散誤差帶來的影響，計(jì)算成本依然過高。此外，基于有限元理論的不確定性參數(shù)建模，均須依托于隨機(jī)概率分布特征和已知概率參數(shù)，導(dǎo)致基于有限元理論的不確定性參數(shù)建模方法使用受限。

小波有限元方法(wavelet finite element method，WFEM)是將傳統(tǒng)有限元方法中使用的多項(xiàng)式插值函數(shù)利用小波函數(shù)代替的有限元方法。由于小波函數(shù)具有多分辨率特征，這就使得在不細(xì)化網(wǎng)格的前提下便可以提高有限元分析方法的求解精度和分辨率，可以同時(shí)兼顧求解精度和求解效率。這使得WFEM自出現(xiàn)以來得到了快速且全面的發(fā)展。依據(jù)文獻(xiàn)資料[7-8]可知，以傳統(tǒng)有限元方法TFEM為代表的確定性分析方法由于計(jì)算成本過高，耗散誤差明顯等原因，存在高頻域難以提供有效數(shù)值解的問題，限制了基于確定性分析方法的寬頻振動(dòng)分析實(shí)現(xiàn)。因此，本文研究聚焦確定性分析方法存在計(jì)算成本過高，耗散誤差明顯等問題，引入具有高求解精度和求解效率特征的WFEM理論開展研究工作，有望突破確定性分析方法無法實(shí)現(xiàn)寬頻振動(dòng)分析的困境。

針對二維問題，Ko等[17]構(gòu)造了基于Daubechies小波的三角形小波單元。Diaz等[18]基于Daubechies小波構(gòu)造了二維中厚板小波單元，并且對該類小波單元的求解能力進(jìn)行了深入研究，研究結(jié)果表明Daubechies小波構(gòu)造的二維中厚板單元不僅可以有效分析中厚板問題，而且可以在對薄板進(jìn)行分析時(shí)避免剪切鎖死現(xiàn)象。為了探究結(jié)合區(qū)間B樣條小波和有限元理論的求解精度和求解效率，Xiang等[19]基于BSWI和單元構(gòu)造理論，構(gòu)造出了適用于求解薄板問題的薄板單元，并在此基礎(chǔ)上展開了多種工況下的數(shù)值分析和實(shí)驗(yàn)研究。而在基于提升框架的小波單元工程應(yīng)用研究領(lǐng)域內(nèi)，Mi等[20]為了提出具有自動(dòng)滿足求精度的有限元分析方法，結(jié)合二代小波函數(shù)理論和有限元理論提出了具有自適應(yīng)網(wǎng)格劃分的二代小波有限元方法，并實(shí)現(xiàn)了對二維聲波問題的求解。對于工程中遇到的奇異性問題求解，Wang等[21]提出了基于提升框架算子的自定義小波有限元方法，并獲得了很好的求解效率和求解精度。在基于Hermite小波的小波有限元方法研究中，Xiang等[22]基于Hermite三次樣條小波構(gòu)造了Hermite小波單元，并在此基礎(chǔ)上對梁結(jié)構(gòu)，轉(zhuǎn)子結(jié)構(gòu)進(jìn)行了模態(tài)分析，相應(yīng)的數(shù)值研究結(jié)果表明該類小波單元具有非常優(yōu)秀的求解精度和求解效率。

可以看出，學(xué)者們基于小波有限元方法開展了豐富的理論和應(yīng)用研究工作，并且針對高模態(tài)密度結(jié)構(gòu)的高頻振動(dòng)分析問題，也開展了部分研究工作，取得了初步的研究成果。本文將簡要介紹依據(jù)小波有限元理論構(gòu)建高模態(tài)密度薄板結(jié)構(gòu)小波有限元分析模型的基本架構(gòu)和寬頻振動(dòng)分析的實(shí)現(xiàn)方式，主要給出了基于小波單元的自耦合算法構(gòu)造架構(gòu)，并給出了基于小波單元和自耦合架構(gòu)實(shí)現(xiàn)高模態(tài)密度結(jié)構(gòu)建模的基本結(jié)構(gòu)，形成了寬頻小波有限元分析方法(wide wavelet finite element method,WWFEM)，著重解決明顯的耗散誤差和計(jì)算成本過高導(dǎo)致傳統(tǒng)有限元方法在對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析時(shí)，難以在高頻域提供精準(zhǔn)的數(shù)值解，致使無法實(shí)現(xiàn)有效的寬頻振動(dòng)分析的問題；隨后結(jié)合數(shù)值分析研究和實(shí)驗(yàn)分析研究方法，討論了小波有限元方法在對高模態(tài)密度薄板結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析時(shí)的有效性和收斂性等，為基于小波有限元方法解決圓柱殼、曲殼等高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析問題提供理論參考。

1 基于薄板理論的小波有限元方法

1.1 小波板單元構(gòu)造

為了基于區(qū)間B樣條小波尺度函數(shù)BSWImj(其中m和j分別為插值小波尺度函數(shù)的階數(shù)和尺度)構(gòu)造如圖2所示的c1BWP單元(小波薄板單元)，首先將如圖2所示的二維求解域Ωe等間隔劃分為n2個(gè)區(qū)域，其中n=2j+m-4。其次將求解域Ωe映射到標(biāo)準(zhǔn)的求解域ΩS，其中Ωs={ξ,η|ξ,η∈[0,1]}。轉(zhuǎn)換到標(biāo)準(zhǔn)求解域Ωs后，各節(jié)點(diǎn)坐標(biāo)可寫為

(2)

圖2 c1BWP單元求解域Fig.2 Solution of c1BWP element

(3)

(4)

式中,Φ1和Φ2均為m階j尺度下的一維BSWI尺度函數(shù)向量。Φ1和Φ2的表達(dá)式如下

(5)

同時(shí)在式(3)中，ae為構(gòu)造二維單元時(shí)的系數(shù)向量，則ae可寫為

(6)

在構(gòu)造c1BWP單元時(shí)，相應(yīng)的物理自由度列向量we可寫為

(7)

隨后，將式(3)代入式(7)可得c1BWP單元的自由度列向量we為

we=Reae

(8)

式(8)中Re可寫為

(9)

(10)

(11)

此時(shí)，式(10)中的形函數(shù)矩陣N可寫為

(12)

因此，式(10)可寫為

(13)

為了構(gòu)造對二維板結(jié)構(gòu)具有自由振動(dòng)分析能力的c1BWP單元，引入板結(jié)構(gòu)自由振動(dòng)勢能泛函Πp，其具體表達(dá)式為

(14)

式中：ρ為板結(jié)構(gòu)材料密度；h為板結(jié)構(gòu)厚度；ω為激振頻率；κ為廣義應(yīng)變矩陣。κ確定表達(dá)式可寫為

(15)

D為板結(jié)構(gòu)的彈性矩陣，可寫為

(16)

式中：D0為板結(jié)構(gòu)彎曲剛度，具體確定表達(dá)式為D0=Eh3/12(1-μ2)，h為板結(jié)構(gòu)厚度，E為材料楊氏模量，μ為泊松比。如式(14)所示，該板結(jié)構(gòu)自由振動(dòng)勢能泛函是在標(biāo)準(zhǔn)求解域Ωs內(nèi)給出的，在此基礎(chǔ)上將式(13)代入式(14)，可得板結(jié)構(gòu)的振動(dòng)模態(tài)方程為

(17)

式中,ω為板結(jié)構(gòu)固有頻率向量。則該類單元對應(yīng)的自由振動(dòng)頻率方程為

(18)

(19)

(20)

1.2 構(gòu)造小波板單元耦合算法

根據(jù)有限元分析理論可知，在建立如圖3(b)所示板子結(jié)構(gòu)小波有限元模型之間的耦合關(guān)系時(shí)，子結(jié)構(gòu)小波有限元模型之間需同時(shí)滿足轉(zhuǎn)角自由度連續(xù)和位移自由度連續(xù)條件。而ID(Index Destination)矩陣主要用于表示全局節(jié)點(diǎn)編號(hào)及其自由度之間的相對關(guān)系，基于該相對關(guān)系可以建立模擬單元之間連續(xù)條件的轉(zhuǎn)換矩陣，從而滿足單元之間的位移自由度連續(xù)條件和轉(zhuǎn)角自由度連續(xù)條件[24]。據(jù)此，首先需建立自耦合索引矩陣(self-coupling index destination，SID)，該矩陣主要用于建立如圖3(b)所示板子結(jié)構(gòu)小波有限元模型在聯(lián)結(jié)處節(jié)點(diǎn)自由度之間的耦合關(guān)系，使得相同節(jié)點(diǎn)處的自由度相等，從而保證板子結(jié)構(gòu)小波有限元模型之間滿足轉(zhuǎn)角自由度和位移自由度連續(xù)性條件。在得到有效SID矩陣的基礎(chǔ)上可建立如圖3(a)所示板結(jié)構(gòu)的整體小波有限元分析模型。

(a) 高模態(tài)密度板結(jié)構(gòu)

(21)

式中：符號(hào)[·]表示對變量x進(jìn)行四舍五入計(jì)算。分段函數(shù)f(x)可基于局部節(jié)點(diǎn)編號(hào)變量lnnode和單元編號(hào)變量ne獲取全局節(jié)點(diǎn)編號(hào)變量gnnode的值。而分段函數(shù)g(x)可基于全局節(jié)點(diǎn)編號(hào)gnnode計(jì)算得到全局自由度編號(hào)gndof。需要聲明的是，在推導(dǎo)SID矩陣過程中同時(shí)還需給出自耦合索引單元節(jié)點(diǎn)矩陣(self-coupling index element nodes，SIEN)，用以表示單元編號(hào)為ne的單元中，局部節(jié)點(diǎn)編號(hào)為lnnode的節(jié)點(diǎn)與全局節(jié)點(diǎn)編號(hào)gnnode之間的對應(yīng)關(guān)系，從而便于建立c1BWP單元之間的對應(yīng)關(guān)系，即全局節(jié)點(diǎn)編號(hào)gnnode與局部節(jié)點(diǎn)編號(hào)lnnode之間的對應(yīng)關(guān)系。

因此，為了建立全局自由度編號(hào)gndof和全局節(jié)點(diǎn)編號(hào)gnnode之間的對應(yīng)關(guān)系，使得子結(jié)構(gòu)小波有限元模型之間滿足位移和轉(zhuǎn)角連續(xù)性條件，必須給出全局節(jié)點(diǎn)編號(hào)gnnode與局部節(jié)點(diǎn)編號(hào)lnnode之間的變化關(guān)系。為了建立該變化關(guān)系的解析表達(dá)式，需構(gòu)造函數(shù)DISdof(a)，該函數(shù)可基于局部節(jié)點(diǎn)編號(hào)lnnode給出相應(yīng)的自由度數(shù)ndof?？梢钥闯?，包括c1BWP單元在內(nèi)的小波板單元根據(jù)局部節(jié)點(diǎn)對應(yīng)的自由度數(shù)均可分為三類。據(jù)此可構(gòu)造具有相同自由度數(shù)的節(jié)點(diǎn)集合DOFs,k，其中k用以表示該集合內(nèi)節(jié)點(diǎn)自由度數(shù)。集合DOFs,k主要用于表示局部節(jié)點(diǎn)編號(hào)lnnode對應(yīng)的自由度數(shù)。因此，c1BWP單元對應(yīng)的集合DOFs,k可寫為

DOFs,4={1，n，n2-n+1，n2}

DOFs,2={a1}∪{a2}∪{a3}∪{a4}

(22)

(23)

與此同時(shí)，式中向量a1，a2，a3和a4可寫為

a1={N∈N+|N=1+i,1≤i≤n-2,i∈N+}

a2={N∈N+|N=1+in,1≤i≤n-2,i∈N+}

a3={N∈N+|N=in,1≤i≤n-2,i∈N+}

a4={N∈N+|N=n2-n+i+1,1≤i≤

n-2,i∈N+}

(24)

此時(shí)，依據(jù)式(22)～式(24)，函數(shù)DISdof(lnnode)可寫為

(25)

可以看出分段函數(shù)DISdof(lnnode)可用以表征局部節(jié)點(diǎn)編號(hào)lnnode與相應(yīng)自由度數(shù)ndof之間的關(guān)系。由上文SID矩陣特征可知，為了建立全局節(jié)點(diǎn)編號(hào)gnnode及其自由度數(shù)之間的對于關(guān)系，必須基于表達(dá)式給出局部節(jié)點(diǎn)編號(hào)lnnode與全局節(jié)點(diǎn)編號(hào)gnnode之間的對于關(guān)系。為此定義函數(shù)R(gnnode)，且lnnode=R(gnnode)，其中R(gnnode)具體形式如下

(26)

式中：n為構(gòu)造小波板單元尺度函數(shù)特征參數(shù)；中間變量A和C確定的表達(dá)式可寫為

(nex-ex+1)(n-1)

(27)

(28)

式中：ex為在水平方向上c1BWP單元個(gè)數(shù)，可以明顯看出，基于函數(shù)R(gnnode)可以建立局部節(jié)點(diǎn)編號(hào)lnnode和全局節(jié)點(diǎn)編號(hào)gnnode之間的內(nèi)在聯(lián)系?；谏鲜鐾茖?dǎo)結(jié)果可以直接建立有效的SID矩陣，與ID矩陣類似，SID矩陣及其元素可以寫為

(29)

此時(shí)，元素sidi,j的求解表達(dá)式可寫為

(30)

可以看出，該式可以通過解析表達(dá)式給出矩陣SID的所有元素，從而建立SID矩陣。除此之外，如上文中所述SIEN主要用于表征編號(hào)為ne的小波單元局部節(jié)點(diǎn)編號(hào)為lnnode的節(jié)點(diǎn)與全局節(jié)編號(hào)gnnode之間的關(guān)系，基于此，SIEN矩陣元素SIEN(lnnode,nelement)可寫為

(31)

(32)

式中,元素si,j可寫為

(33)

SIEN(lnnode,ne)}

(34)

(35)

式中：n表征了建立高模態(tài)密度結(jié)構(gòu)引用的小波單元類型，本章引用的是c1BWP單元，相應(yīng)的插值小波尺度函數(shù)階數(shù)和尺度分為4和3，此時(shí)n為9。使用其他小波單元建立高模態(tài)密度結(jié)構(gòu)小波有限元分析模型時(shí)可依據(jù)實(shí)際情況進(jìn)行修正。此時(shí)，基于如圖3(b)所示子結(jié)構(gòu)劃分方式，編號(hào)為ne的小波單元c1BWP與板結(jié)構(gòu)小波有限元數(shù)值模型之間的耦合關(guān)系可寫為

(36)

(37)

(2) 將板結(jié)構(gòu)小波有限元模型的整體剛度矩陣Khigh和整體質(zhì)量矩陣Μhigh代入至模態(tài)方程

(Khigh-ω2Mhigh)we=0

(38)

并建立相應(yīng)的特征方程

|Khigh-ω2Mhigh|=0

(39)

從而計(jì)算得到相應(yīng)的頻率向量ω和模態(tài)振型Φ，其中ω和Φ可寫為

Φ=[φ1,φ2,…,φi,…,φn]

ω=[ω1,ω2,…,ωi,…ωn]

(40)

(41)

(42)

圖4 寬頻小波有限元分析方法(WWFEM)開展高模態(tài)密度板結(jié)構(gòu)寬頻振動(dòng)分析流程圖Fig.4 Chart of WWFEM for analyzing the high modal plate

上述即為基于寬頻小波有限元方法對高模態(tài)密度板結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的基本過程，隨后，將分別基于數(shù)值分析研究和實(shí)驗(yàn)分析研究，討論該分析方法的有效性。

2 寬頻振動(dòng)的小波有限元方法數(shù)值分析

2.1 有效性分析

為了基于量化指標(biāo)定義高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析，從而便于開展高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析方法驗(yàn)證研究。采用我國著名學(xué)者姚德源在其專著[25]中給出的依據(jù)帶寬內(nèi)的模態(tài)疊合因子(modal overlap factor，MOF)，即帶寬內(nèi)的模態(tài)數(shù)，來作為量化指標(biāo)的寬頻振動(dòng)分析定義方法，其將寬頻域劃分低頻域，中頻域和高頻域，具體如下

(43)

并將在上述頻域內(nèi)進(jìn)行的振動(dòng)分析分別定義為低頻振動(dòng)分析，中頻振動(dòng)分析和高頻振動(dòng)分析，這也是本文進(jìn)行高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析方法驗(yàn)證研究的基本定義。

如前文所述，當(dāng)確定性分析方法存在計(jì)算成本過高，耗散誤差明顯時(shí)，將導(dǎo)致其在對高模態(tài)密度板結(jié)構(gòu)進(jìn)行振動(dòng)分析時(shí)難以在較高頻域內(nèi)提供有效數(shù)值解，導(dǎo)致無法對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析。因此，如果WWFEM可在高頻域內(nèi)可對高模態(tài)密度結(jié)構(gòu)實(shí)現(xiàn)有效的振動(dòng)分析，則可認(rèn)為WWFEM必然具有可同時(shí)在低頻域，中頻域和高頻域?qū)Ω吣B(tài)密度結(jié)構(gòu)進(jìn)行振動(dòng)分析的求解能力，隨后將基于該思路驗(yàn)證小波有限元方法對高模態(tài)密度結(jié)構(gòu)的寬頻振動(dòng)分析能力。

為此，引入如圖5所示的四邊簡支薄板結(jié)構(gòu)數(shù)值模型，其中l(wèi)x和ly分別為該結(jié)構(gòu)在水平和垂直方向的長度，該數(shù)值模型的幾何，物理參數(shù)如表1所示?；诖?，首先可以看出該板結(jié)構(gòu)數(shù)值模型厚度h與特征尺寸(矩形板對角線長度)的比值遠(yuǎn)遠(yuǎn)小于1/10，依據(jù)高模態(tài)密度結(jié)構(gòu)定義可知其為典型的高模態(tài)密度結(jié)構(gòu)。因此，該數(shù)值模型可用于分析研究WWFEM對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的有效性?；趯Ρ确治龇椒?yàn)證WWFEM有效性的參考解是基于前1 600階模態(tài)振型和固有頻率解析解得到的振動(dòng)響應(yīng)解析解。

圖5 四邊簡支高模態(tài)密度板結(jié)構(gòu)數(shù)值模型示意圖Fig.5 Diagram of numerical model of high-mode density plate with simply supported on four sides

表1 高模態(tài)密度板結(jié)構(gòu)材料參數(shù)Tab.1 Parameters of high modal density plate

依據(jù)如式(43)所示的高模態(tài)密度結(jié)構(gòu)寬頻域劃分方式，結(jié)合1/3倍頻程將分析頻域Ω為5～1 000 Hz的頻域劃分為若干個(gè)分析帶寬，在此基礎(chǔ)上結(jié)合如圖5所示數(shù)值模型的固有頻率解析解計(jì)算各個(gè)帶寬內(nèi)的模態(tài)數(shù)MOF?？傻迷摂?shù)值模型帶寬內(nèi)的模態(tài)數(shù)隨著帶寬編號(hào)的分布特征，具體如圖6所示，其中(1,0)和(23,45)中，括號(hào)左邊代表帶寬編號(hào)，右邊代表帶寬內(nèi)的模態(tài)數(shù)。由圖中可以看出，在分析頻域內(nèi)該高模態(tài)密度結(jié)構(gòu)數(shù)值模型在帶寬內(nèi)的模態(tài)數(shù)最高為45，即MOF最高為45?？梢?，該數(shù)值模型在分析頻域內(nèi)的MOF遠(yuǎn)遠(yuǎn)大于式(43)所示的高頻域閾值5。因此，如果小波有限元方法可在分析頻域?yàn)?～1 000 Hz內(nèi)對如圖5所示數(shù)值模型進(jìn)行有效的振動(dòng)分析，則依據(jù)上述寬頻振動(dòng)分析能力驗(yàn)證思路，可認(rèn)為WWFEM具有對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的能力。

圖6 帶寬內(nèi)的模態(tài)數(shù)MOF隨帶寬分布Fig.6 Distribution of MOF with bandwidth

據(jù)此，分別基于WWFEM和解析方法求解如圖5所述的數(shù)值模型在5～1 000 Hz以內(nèi)的加速度響應(yīng)，其中激振點(diǎn)位置為(xe,ye)=(0.125 m,0.25 m)，拾振點(diǎn)位置為(xr,yr)=(1.375 m,0.75 m)。得到的加速度數(shù)值解和解析解如圖7所示，從圖中可以看出，WWFEM在分析頻域Ω為5～1 000 Hz內(nèi)提供的加速度數(shù)值解與解析解保持非常高的一致性。同時(shí)依據(jù)上文所述寬頻振動(dòng)分析能力驗(yàn)證思路可見，WWFEM可以有效地對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析。

(a)

為了進(jìn)一步驗(yàn)證WWFEM對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的有效性，隨后將研究分析WWFEM解決現(xiàn)有確定性方法存在的計(jì)算成本過高和耗散誤差明顯問題的能力，為此將分別討論小波有限元方法的求解精度(可以用以反映耗散誤差的大小)和求解效率。

2.2 WWFEM求解精度和效率研究

依據(jù)WWFEM方法寬頻振動(dòng)分析研究結(jié)果可以看出，WWFEM方法可對高模態(tài)密度結(jié)構(gòu)進(jìn)行有效地寬頻振動(dòng)分析。本部分將主要對其能夠突破現(xiàn)有確定性分析方法局限性，基于數(shù)值分析方法對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的主要原因進(jìn)行討論，即討論分析其解決現(xiàn)有確定性方法存在的計(jì)算成本過高和耗散誤差明顯問題的能力。為此，討論主要分為兩方面，一方面討論WWFEM的求解精度是否可以解決現(xiàn)有確定性分析方法存在耗散誤差明顯的問題，另一方面討論WWFEM的求解效率是否可以解決現(xiàn)有確定性分析方法計(jì)算成本過高的問題。

為了說明WWFEM可以解決現(xiàn)有確定性分析方法存在耗散誤差明顯的問題，仍以圖5所示高模態(tài)密度板結(jié)構(gòu)為數(shù)值模型對WWFEM的求解精度進(jìn)行對比分析，相關(guān)參數(shù)同上。

為了實(shí)現(xiàn)對比，將傳統(tǒng)有限元方法(traditional finite element method,TFEM)方法作為確定性分析方法參考解，并基于Ansys作為TFEM方法的實(shí)現(xiàn)平臺(tái)提供確定性分析方法參考數(shù)值解。據(jù)此，首先利用Ansys平臺(tái)對圖5所示數(shù)值模型進(jìn)行建模，隨后在不同網(wǎng)格劃分條件下獲取激振點(diǎn)位置為(xe,ye)=(0.125 m,0.25 m)，拾振點(diǎn)位置為(xr,yr)=(1.375 m,0.75 m)的加速度響應(yīng)解，從而獲取Ansys最高求解精度條件下對應(yīng)的振動(dòng)特征數(shù)值解。研究發(fā)現(xiàn)，當(dāng)網(wǎng)格劃分條件分別為8×12,80×120和128×192時(shí)，對應(yīng)在200～300 Hz內(nèi)的穩(wěn)態(tài)響應(yīng)解如圖8所示。

圖8 Ansys(Shell63)求解精度極限Fig.8 Solution accuracy of Ansys(Shell63)

從圖中可以看出，當(dāng)網(wǎng)格數(shù)為128×192時(shí)Ansys計(jì)算平臺(tái)的計(jì)算精度已經(jīng)難以再通過網(wǎng)格細(xì)化進(jìn)行有效提升。因此可見，當(dāng)網(wǎng)格劃分條件為128×192時(shí)300 Hz以內(nèi)的振動(dòng)分析數(shù)值解即為Ansys能夠給出的最高求解精度。

為了基于量化對比分析方法研究分析WWFEM與Ansys最高求解精度的差異，以300 Hz以內(nèi)的前57階固有頻率展開定量對比分析研究。為了可以直觀分析WWFEM和Ansys求解精度的差異，定義了相對誤差ε(i)，計(jì)算表達(dá)式為

(44)

式中：ωN(i)為基于WWFEM方法和Ansys得到的第i階固有頻率數(shù)值解；ωA(i)為基于解析表達(dá)式得到了第i階固有頻率解析解。將基于WWFEM方法和Ansys得到的固有頻率數(shù)值解代入式(44)可得前57階固有頻率相對誤差，具體結(jié)果如圖9所示，圖中實(shí)線表示為Ansys在預(yù)測前57解固有頻率時(shí)的最高精度，虛線為WWFEM方法在預(yù)測前57階固有頻率時(shí)得到的數(shù)值解對應(yīng)的相對誤差?？梢悦黠@看出，Ansys最高精度的誤差均明顯高于WWFEM的求解誤差，且求解時(shí)間約為。除此之外，如圖9所示并不是WWFEM方法的求解精度極限，該求解精度仍然可再次提升，甚至隨著劃分板子結(jié)構(gòu)數(shù)量的提升，數(shù)值解可與解析解保持完全一致，而TFEMs由于耗散誤差明顯難以做到這一點(diǎn)。同時(shí)觀察Ansys給出的最高求解精度可以看出，隨著模態(tài)數(shù)的增加，其求解精度呈現(xiàn)明顯降低且不斷波動(dòng)的趨勢，此外結(jié)合圖8所示可以看出，上述精度已經(jīng)是Ansys方法的計(jì)算極限，其求解精度已無法再次提升，而隨著頻域范圍越來越高，分析帶寬內(nèi)需要精準(zhǔn)計(jì)算的模態(tài)數(shù)將急劇增加，這使得Ansys已經(jīng)無法實(shí)現(xiàn)更高頻域的振動(dòng)特征分析，究其原因主要是由于TFEMs明顯的耗散誤差。而WWFEM成功地突破了TFEMs理論的求解精度限制，隨著模態(tài)數(shù)增加求解精度變化穩(wěn)定，誤差退化較慢，此外該誤差仍然可以通過劃分更多的板子結(jié)構(gòu)模型實(shí)現(xiàn)穩(wěn)定降低。并且WWFEM在獲取上述數(shù)值結(jié)果過程中的計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)小于TFEMs，約為15 s左右。

圖9 WWFEM求解精度對比Fig.9 Comparison of WWFEM solution accuracy

此外，如表2所示為WWFEM和TFEM在獲得同樣求解精度時(shí)相應(yīng)的建模和計(jì)算成本，其中Nnodes為模型節(jié)點(diǎn)數(shù)，Ndofs為模型自由度數(shù)?？梢钥闯霎?dāng)基于TFEM方法得到前100階固有頻率，數(shù)值解最大相對誤差為0.3%時(shí)，相應(yīng)的高模態(tài)密度板結(jié)構(gòu)有限元模型節(jié)點(diǎn)數(shù)為14 095，相應(yīng)的自由度數(shù)為84 570左右。而WWFEM得到前100階固有頻率數(shù)值解最大相對誤差為0.3%時(shí)，高模態(tài)密度板結(jié)構(gòu)小波有限元模型節(jié)點(diǎn)數(shù)僅為1 089，而自由度僅為1 316?？梢钥闯鯳WFEM建模包含的自由度數(shù)僅為傳統(tǒng)有限元的1/64，這從側(cè)面可以反映出WWFEM的計(jì)算成本遠(yuǎn)低于TFEM。除此之外，WWFEM在達(dá)到求解精度為0.3%時(shí)的計(jì)算時(shí)間小于4 s，而TFEMs方法需要的計(jì)算時(shí)間將明顯大于10 000 s(約3 h)。

表2 WWFEM方法和TFEM求精效率對比分析Tab.2 Comparative analysis of WWFEM and TFEM

可見，WWFEM可以解決TFEM存在的耗散誤差明顯和成本過高的問題，從而突破TFEM的求解能力極限，補(bǔ)足了確定性分析方法對高模態(tài)密度結(jié)構(gòu)進(jìn)行高頻振動(dòng)分析時(shí)難以滿足精度要求的缺陷。因此，WWFEM可對高模態(tài)密度結(jié)構(gòu)實(shí)現(xiàn)有效的寬頻振動(dòng)分析。為了更進(jìn)一步研究WWFEM方法的有效性，隨后將開展實(shí)驗(yàn)分析研究。

3 寬頻振動(dòng)的小波有限元方法實(shí)驗(yàn)分析

3.1 實(shí)驗(yàn)研究

為了基于實(shí)驗(yàn)分析研究WWFEM對高模態(tài)密度板結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的有效性，搭建了如圖10所示的板結(jié)構(gòu)實(shí)驗(yàn)平臺(tái)，相應(yīng)的示意圖如圖11所示。從該示意圖中可以看出，實(shí)驗(yàn)研究的板結(jié)構(gòu)邊界條件為兩長邊固支，兩短邊自由，相應(yīng)的幾何，物理參數(shù)如表3所示。由此可知該板結(jié)構(gòu)特征尺寸(矩形板對角線長度)與厚度比值遠(yuǎn)遠(yuǎn)小于1/10，即該結(jié)構(gòu)為典型的高模態(tài)密度結(jié)構(gòu)。因此，基于圖10所示板結(jié)構(gòu)開展的振動(dòng)實(shí)驗(yàn)研究可用于分析WWFEM方法對高模態(tài)密度結(jié)構(gòu)的寬頻振動(dòng)分析能力。

(a)

圖11 高模態(tài)密度板實(shí)驗(yàn)?zāi)Ｐ褪疽鈭DFig.11 Schematic diagram of high mode density plate

表3 高模態(tài)密度板結(jié)構(gòu)實(shí)物幾何材料參數(shù)Tab.3 Parameters of high modal density board structure

而為了基于實(shí)驗(yàn)研究方法分析WWFEM方法對高模態(tài)密度板結(jié)構(gòu)的寬頻振動(dòng)分析能力，必須在寬頻域內(nèi)對比分析WWFEM提供的數(shù)值解與實(shí)驗(yàn)解的相似性，從而說明其具有寬頻振動(dòng)分析的能力。

為此，需確定寬頻域Ω，使其同時(shí)包含有如式(43)所示的低頻域，中頻域和高頻域。為此，首先將頻域Ω為5～1 000 Hz基于1/3倍頻程劃分為若干個(gè)帶寬Δωi，進(jìn)而基于式(40)得到的固有頻率向量得到不同帶寬內(nèi)的模態(tài)數(shù)MOF，得到不同帶寬內(nèi)的模態(tài)數(shù)分布特征如圖12所示。從該圖中可以明顯看出，不同帶寬Δωi內(nèi)的模態(tài)數(shù)在0～22之間變化，即如圖10所示的高模態(tài)密度板結(jié)構(gòu)實(shí)驗(yàn)平臺(tái)在5～1 000 Hz頻域內(nèi)MOF最大為22，明顯大于如式(43)所示的高頻域閾值5，因此Ω為5～1 000 Hz同時(shí)包含有低頻域，中頻域和高頻域。換言之，在頻域Ω為5～1 000 Hz時(shí)對圖10所示高模態(tài)密度板結(jié)構(gòu)實(shí)驗(yàn)平臺(tái)進(jìn)行的振動(dòng)分析時(shí)，同時(shí)包含有低頻振動(dòng)分析，中頻振動(dòng)分析和高頻振動(dòng)分析。因此，對如圖10所示高模態(tài)密度板結(jié)構(gòu)在頻域Ω內(nèi)進(jìn)行的振動(dòng)分析實(shí)驗(yàn)研究，可以研究分析WWFEM的寬頻振動(dòng)分析能力。

圖12 模態(tài)數(shù)分布特征Fig.12 Distribution characteristics of modal number

據(jù)此，首先分別基于實(shí)驗(yàn)方法和WWFEM得到如圖10所示高模態(tài)密度板結(jié)構(gòu)在頻域Ω內(nèi)的加速度響應(yīng)數(shù)值解和實(shí)驗(yàn)結(jié)果，隨后對比得到的加速度響應(yīng)數(shù)值解和實(shí)驗(yàn)結(jié)果，依據(jù)對比結(jié)果研究分析WWFEM對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的有效性。

為了實(shí)現(xiàn)上述目標(biāo)，實(shí)驗(yàn)過程采用的儀器主要包括：力錘，加速度傳感器和信號(hào)采集系統(tǒng)，儀器型號(hào)和實(shí)驗(yàn)時(shí)選擇的靈敏度信息如表4所示。與此同時(shí)，為了滿足數(shù)據(jù)精度需要，在實(shí)驗(yàn)過程中設(shè)置的采樣率為2 560 Hz，對應(yīng)的分辨率為0.071 85 Hz。實(shí)驗(yàn)時(shí)，將傳感器依據(jù)表5所示的拾振點(diǎn)位置信息安裝于如圖10所示的高模態(tài)密度板結(jié)構(gòu)振動(dòng)實(shí)驗(yàn)臺(tái)中，隨后依據(jù)表5所示的激振點(diǎn)位置信息使用力錘在相應(yīng)位置點(diǎn)處進(jìn)行敲擊實(shí)驗(yàn)。最后，將采集到的加速度信號(hào)和力信號(hào)輸入至信號(hào)處理系統(tǒng)可得加速度響應(yīng)實(shí)驗(yàn)結(jié)果。

表4 實(shí)驗(yàn)設(shè)備信息Tab.4 Laboratory equipment information

表5 激振點(diǎn)和拾振點(diǎn)位置信息Tab.5 Information of excitation point and pick-up point

3.2 WWFEM振動(dòng)實(shí)驗(yàn)結(jié)果分析

依據(jù)3.1中所述可知，為了基于實(shí)驗(yàn)對比分析方法研究分析WWFEM對高模態(tài)密度結(jié)構(gòu)的寬頻振動(dòng)分析能力。首先，基于錘擊法和表5所示激振點(diǎn)位置信息和拾振位置信息得到加速度響應(yīng)的實(shí)驗(yàn)結(jié)果。隨后，基于WWFEM和如圖11所示板結(jié)構(gòu)實(shí)驗(yàn)平臺(tái)的幾何，材料參數(shù)(如表3所示)建立如圖10所示高模態(tài)密度板結(jié)構(gòu)的小波有限元分析模型，并在此基礎(chǔ)上依據(jù)表5所示激振點(diǎn)位置信息和拾振位置信息得到加速度響應(yīng)數(shù)值解。最后，依據(jù)得到的加速度響應(yīng)數(shù)值解和實(shí)驗(yàn)結(jié)果展開對比分析研究，具體對比分析結(jié)果如圖13所示。

結(jié)合3.1中寬頻域定義可以看出，在5～300 Hz之間的較低頻域內(nèi)，自左往右觀察圖13(a)～圖13(d)所示數(shù)值解與實(shí)驗(yàn)解可以看出兩者均保持非常好的一致性(在聲振分析研究領(lǐng)域)，這說明了WWFEM方法在較低頻域內(nèi)對高模態(tài)密度進(jìn)行振動(dòng)分析的可靠性。

與此同時(shí)，自左往右觀察如圖13 (a)和圖13 (c)所示結(jié)果可以看出，數(shù)值解在頻域?yàn)?00～1 000 Hz之間仍然與實(shí)驗(yàn)解保持著高度的一致性(在聲振分析研究領(lǐng)域)，而該頻域明顯包含有高頻域的。因此WWFEM具有自低頻到高頻的振動(dòng)分析能力。依據(jù)式(43)所述可見，WWFEM方法具有可對高模態(tài)密度板結(jié)構(gòu)進(jìn)行由低頻到高頻的寬頻振動(dòng)分析能力。而圖13(b)和圖13(d)所示數(shù)值解在頻域?yàn)?00～1 000 Hz之間與實(shí)驗(yàn)解的誤差明顯高于圖13(a)和圖13(c)所示，該誤差主要是由于在對應(yīng)的實(shí)驗(yàn)中錘擊法在較高頻域內(nèi)較低的信噪比造成的。除此之外，WWFEM在得到如圖13所示的高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)響應(yīng)數(shù)值解時(shí)基于個(gè)人PC計(jì)算平臺(tái)的計(jì)算時(shí)間均小于15 s。因此可見，WWFEM對實(shí)際工程中的高模態(tài)密度板結(jié)構(gòu)可以實(shí)現(xiàn)有效的寬頻振動(dòng)分析，兼顧高求解效率。

(a) 工況1

此外，為了進(jìn)一步說明本文所提方法進(jìn)行寬頻振動(dòng)分析時(shí)的先進(jìn)性。引入Hybrid SEA/TFEMs架構(gòu)中對高模態(tài)密度結(jié)構(gòu)實(shí)現(xiàn)振動(dòng)分析的SEA方法和本文提出的WWFEM方法開展寬頻振動(dòng)分析驗(yàn)證研究，得到的對比分析結(jié)果如圖14所示。自左往右觀察如圖14所示結(jié)果可以看出，在20～2 000 Hz的頻域范圍內(nèi)，由WWFEM計(jì)算得到的加速度響應(yīng)數(shù)值解與實(shí)驗(yàn)得到的加速度響應(yīng)實(shí)驗(yàn)結(jié)果保持著非常好的一致性，尤其在高頻域1 500 Hz附近提供了精準(zhǔn)的振動(dòng)相應(yīng)細(xì)節(jié)特征。而Hybrid SEA/TFEMs架構(gòu)下的SEA方法雖然在高頻域可以提供振動(dòng)水平的平均情況，但仍然存在低頻域誤差較大，并且在高頻域存在無法提供振動(dòng)響應(yīng)細(xì)節(jié)信息的問題。

圖14 加速度頻響函數(shù)數(shù)值解對比分析Fig.14 Comparison of experimental numerical solutions

綜上所述可見WWFEM可以僅基于確定性分析方法提供比SEA方法更為有效的數(shù)值解，同時(shí)非常有效地解決了SEA方法在進(jìn)行寬頻振動(dòng)分析時(shí)存在的低頻域誤差較大、無法提供振動(dòng)響應(yīng)細(xì)節(jié)信息的問題。這為簡化振動(dòng)分析方法的使用提供了有效的技術(shù)支撐，有望僅基于有限元分析平臺(tái)即可實(shí)現(xiàn)高模態(tài)密度結(jié)構(gòu)的寬頻域振動(dòng)分析。

4 結(jié) 論

以傳統(tǒng)有限元分析方法為代表的確定性分析方法由于計(jì)算成本過高，存在明顯耗散誤差等難題，使得基于傳統(tǒng)有限元理論的商業(yè)軟件存在難以對高模態(tài)密度結(jié)構(gòu)進(jìn)行寬頻振動(dòng)分析的問題。本文結(jié)合理論推導(dǎo)、數(shù)值分析驗(yàn)證和實(shí)驗(yàn)分析驗(yàn)證，介紹了小波有限元方法在解決該類問題時(shí)的潛在優(yōu)勢。并重點(diǎn)論述了自耦合算法的推導(dǎo)架構(gòu)，為基于小波單元建立高模態(tài)密度結(jié)構(gòu)分析模型提供理論基礎(chǔ)，形成了寬頻小波有限元分析方法，并對該方法的有效性進(jìn)行了數(shù)值分析研究和實(shí)驗(yàn)分析研究。結(jié)果顯示，在建模精度符合要求的條件下，小波有限元分析方法可對高模態(tài)密度板結(jié)構(gòu)進(jìn)行非常有效的寬頻振動(dòng)分析，并且可快速提供數(shù)值解，保持優(yōu)秀的穩(wěn)定性。本文研究結(jié)果可為基于小波有限元分析方法解決圓柱殼和曲殼等經(jīng)典高模態(tài)密度結(jié)構(gòu)的寬頻振動(dòng)分析問題提供理論依據(jù)，對具有復(fù)雜幾何特征的高模態(tài)密度結(jié)構(gòu)寬頻振動(dòng)分析，仍需進(jìn)一步研究。