基于功率流有限元方法的異形薄板能量密度求解

劉知輝， 牛軍川,2， 周一群

(1.山東大學(xué) 機(jī)械工程學(xué)院，濟(jì)南 250061；2. 山東大學(xué) 高效潔凈機(jī)械制造教育部重點(diǎn)實(shí)驗(yàn)室，濟(jì)南 250061)

基于功率流有限元方法的異形薄板能量密度求解

劉知輝1， 牛軍川1,2， 周一群1

(1.山東大學(xué) 機(jī)械工程學(xué)院，濟(jì)南 250061；2. 山東大學(xué) 高效潔凈機(jī)械制造教育部重點(diǎn)實(shí)驗(yàn)室，濟(jì)南 250061)

在機(jī)械工程中，對(duì)結(jié)構(gòu)在不同頻率激勵(lì)下的振動(dòng)響應(yīng)進(jìn)行分析預(yù)測(cè)具有重要的意義。功率流有限元法以其適用頻率范圍較廣，可得到結(jié)構(gòu)響應(yīng)的細(xì)節(jié)信息等優(yōu)點(diǎn)成為振動(dòng)分析領(lǐng)域的研究熱點(diǎn)。利用功率流有限元方法對(duì)薄板的彎曲波能量密度進(jìn)行了求解，使用加權(quán)殘差法導(dǎo)出了薄板單元節(jié)點(diǎn)的能量密度殘差，利用線性四邊形網(wǎng)格對(duì)薄板進(jìn)行網(wǎng)格劃分并在此基礎(chǔ)上建立了單元的有限元方程，進(jìn)一步地通過(guò)對(duì)單元有限元方程的組裝和求解得到了薄板上各個(gè)節(jié)點(diǎn)處的能量密度響應(yīng)，引入線性三角形網(wǎng)格以處理復(fù)雜形狀薄板能量密度的求解，對(duì)使用功率流有限元方法求解任意形狀薄板上的能量密度分布問(wèn)題具有一定意義。

振動(dòng)；功率流有限元分析；異形薄板；能量密度

結(jié)構(gòu)的振動(dòng)對(duì)結(jié)構(gòu)的性能、壽命有重要影響，對(duì)于給定的結(jié)構(gòu)，在設(shè)計(jì)階段使用有效的分析方法對(duì)結(jié)構(gòu)的振動(dòng)響應(yīng)進(jìn)行分析并得到較為準(zhǔn)確的結(jié)果對(duì)于結(jié)構(gòu)的設(shè)計(jì)和應(yīng)用具有重要的意義[1]。

傳統(tǒng)的結(jié)構(gòu)振動(dòng)分析方法主要有基于位移的有限元方法和基于能量的統(tǒng)計(jì)能量法，其中有限元方法在低頻時(shí)可以較為準(zhǔn)確地預(yù)測(cè)結(jié)構(gòu)的振動(dòng)響應(yīng)，但當(dāng)頻率增大時(shí)，為了準(zhǔn)確得到結(jié)構(gòu)的響應(yīng)不得不對(duì)網(wǎng)格進(jìn)行細(xì)化，由此帶來(lái)的不斷增加的運(yùn)算量最終將使得計(jì)算變得不可行，且在高頻時(shí)結(jié)構(gòu)的響應(yīng)對(duì)結(jié)構(gòu)的細(xì)節(jié)變得極為敏感，表現(xiàn)出混沌的現(xiàn)象，有限元建模過(guò)程中對(duì)細(xì)節(jié)的忽略以及數(shù)據(jù)的不確定性使得在高頻時(shí)分析得到的結(jié)果不具有可信性，因此，在高頻時(shí)使用有限元方法難以有效的預(yù)測(cè)結(jié)構(gòu)的響應(yīng)[2]。統(tǒng)計(jì)能量法在分析結(jié)構(gòu)的高頻振動(dòng)時(shí)具有較大的優(yōu)勢(shì)，可以較為準(zhǔn)確的預(yù)測(cè)高頻時(shí)結(jié)構(gòu)在統(tǒng)計(jì)意義下的響應(yīng)。然而由于統(tǒng)計(jì)能量法在求解結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)時(shí)一般將具體的耦合結(jié)構(gòu)劃分為若干個(gè)子系統(tǒng)，最終求解得到每個(gè)子系統(tǒng)在統(tǒng)計(jì)意義下的平均能量響應(yīng)，無(wú)法再提供子系統(tǒng)內(nèi)部細(xì)節(jié)信息，因此當(dāng)需要了解子系統(tǒng)內(nèi)局部的響應(yīng)時(shí)，統(tǒng)計(jì)能量法便不再適用[3]。由Nefske等[4]提出的功率流有限元方法將結(jié)構(gòu)中每一微元視為一子系統(tǒng)，利用能量平衡關(guān)系建立結(jié)構(gòu)的能量密度控制方程，并進(jìn)一步的通過(guò)求解能量密度控制方程得到結(jié)構(gòu)的能量密度響應(yīng)。與傳統(tǒng)的有限元方法相比，該方法可較為有效的分析結(jié)構(gòu)的中高頻振動(dòng)響應(yīng)，與統(tǒng)計(jì)能量法相比，該方法可得到結(jié)構(gòu)內(nèi)局部響應(yīng)的信息，因此該方法在結(jié)構(gòu)中高頻動(dòng)力學(xué)分析領(lǐng)域具有較大的優(yōu)越性。在Nefske等完成了功率流有限元開(kāi)創(chuàng)性的工作后，眾多學(xué)者對(duì)具體結(jié)構(gòu)的能量密度控制方程和耦合結(jié)構(gòu)間的能量傳遞關(guān)系進(jìn)行了研究[5-7]。

板是在工程實(shí)際中常見(jiàn)的結(jié)構(gòu)，在許多耦合系統(tǒng)中板的輻射聲是系統(tǒng)向外輻射聲能的重要組成部分，因此使用功率流有限元方法研究外界激勵(lì)下板的能量密度分布具有重要意義。Bouthier等[8-9]先后研究并建立了薄膜和薄板的彎曲波能量密度控制方程，Park等[10]對(duì)薄板的面內(nèi)波能量密度控制方程進(jìn)行了研究并基于此研究了兩塊矩形板耦合時(shí)能量傳遞特性，李坤朋[11]研究了矩形薄板以及矩形薄板中存在加強(qiáng)筋時(shí)板上能量的分布和傳遞，Niu等[12- 13]研究了考慮三種波形情況下耦合板的能量密度求解，江民圣等[14]研究了兩矩形薄板以任意角度耦合時(shí)板中能量密度的分布與傳遞特性。上述研究對(duì)功率流有限元方法應(yīng)用于結(jié)構(gòu)中板的能量密度求解有著重要意義，但這些研究都以規(guī)則的矩形板為研究對(duì)象，且劃分的單元皆為統(tǒng)一的矩形單元，對(duì)于工程實(shí)際中應(yīng)用更為廣泛的形狀不規(guī)則薄板無(wú)法進(jìn)行仿真計(jì)算，且缺少相關(guān)研究。

本文利用功率流有限元方法對(duì)薄板的能量密度進(jìn)行了求解，利用加權(quán)殘值法建立了單元節(jié)點(diǎn)的能量密度殘差方程，使用線性四邊形單元對(duì)形狀較規(guī)則薄板進(jìn)行網(wǎng)格劃分，求解了薄板的能量密度分布并對(duì)薄板中的能量密度分布和能流進(jìn)行了可視化，進(jìn)一步的利用三角形單元對(duì)復(fù)雜形狀的薄板進(jìn)行網(wǎng)格劃分并建立了三角形單元的能量密度有限元方程，同時(shí)對(duì)不同單元類(lèi)型混合使用時(shí)的單元兼容性進(jìn)行了討論，三角形單元的引入以及不同單元類(lèi)型的混合使用可以有效的將功率流有限元方法應(yīng)用范圍擴(kuò)展至任意幾何形狀的薄板結(jié)構(gòu)。

1 能量密度控制方程

如圖1所示，薄板受到外界的激勵(lì)后產(chǎn)生對(duì)應(yīng)的響應(yīng)，從而獲得動(dòng)能與勢(shì)能，定義能量密度為單位面積內(nèi)動(dòng)能與勢(shì)能的總和，也即

e=T+V

(1)

式中，e、T和V分別為系統(tǒng)的總能量密度、動(dòng)能和勢(shì)能。若薄板處于穩(wěn)態(tài)，則其任意位置處能量密度為常數(shù)，即

(2)

如圖2所示，考慮處于穩(wěn)態(tài)的薄板內(nèi)任意一微元內(nèi)的能量平衡關(guān)系

-∫Γn·qdΓ-pdiss+pin=0

(3)

圖1 受功率激勵(lì)的薄板Fig. 1 Thin plate excited by the power input

式中：q為能流強(qiáng)度；Г為微元的邊界；n為微元單元邊界外法線方向矢量；pdiss為結(jié)構(gòu)阻尼引起的能量耗散；pin為外界的功率輸入。式(3)刻畫(huà)了微元內(nèi)的能量平衡關(guān)系。

圖2 微元內(nèi)的能量平衡Fig. 2 Energy balance in the infinitesimal element

基于平面波遠(yuǎn)場(chǎng)疊加原則和式(3)，Bouthier等分別具體的導(dǎo)出了各向同性有限薄板的彎曲波和面內(nèi)波的能量密度控制方程

(4)

式中：cg為彈性波的群速度；η為結(jié)構(gòu)損耗因子；ω為激勵(lì)圓頻率；Δ(·)二維拉普拉斯算子；e為周期平均以及局部空間平均意義下遠(yuǎn)場(chǎng)波能量密度；πin外界輸入功率。

式(4)確定了能量密度在薄板內(nèi)的分布情況，其中群速度cg由彈性波的具體類(lèi)型確定，當(dāng)彈性波為具有非頻散特性的面內(nèi)彈性波時(shí)，群速度cg與波速c的關(guān)系為cg=c，當(dāng)彈性波為具有頻散特性的板面外彎曲波時(shí)cg= 2c。 Seo等[15]使用單傅里葉級(jí)數(shù)類(lèi)型的Lévy解對(duì)式(4)進(jìn)行求解，該方法對(duì)求解域有較高的要求，難以應(yīng)用于復(fù)雜的求解域。Pereira等[16]使用譜元法對(duì)式(4)進(jìn)行求解，同樣存在不能適用于復(fù)雜求解域的缺陷。有限元可適用于復(fù)雜求解域，并可將結(jié)構(gòu)的局部特性考慮在內(nèi)，因此可考慮使用有限元方法對(duì)式(4)進(jìn)行求解。由于能量密度e的泛函不易獲得，因此難以通過(guò)變分原理建立式(4)有限元方程，可考慮使用以形函數(shù)作為權(quán)函數(shù)的加權(quán)殘差法對(duì)式(4)進(jìn)行求解。

單板結(jié)構(gòu)一般主要考慮其彎曲波的影響，式(4)對(duì)于薄板內(nèi)的彎曲波能量密度控制方程的研究基于薄板的Kirchhoff-Love理論，此時(shí)薄板的變形使用中性層的撓度進(jìn)行描述，且式(4)在建立過(guò)程中已對(duì)厚度方向進(jìn)行積分，因此對(duì)薄板上能量密度的求解類(lèi)似于平面應(yīng)力問(wèn)題的求解，能量密度僅取決于板的面內(nèi)位置而與厚度方向位置無(wú)關(guān)，因此在使用有限元方法對(duì)板上的能量密度進(jìn)行求解時(shí)，應(yīng)使用二維平面網(wǎng)格對(duì)薄板的幾何區(qū)域進(jìn)行離散。在形函數(shù)加權(quán)下，離散后所得到的二維平面網(wǎng)格單元Ω的節(jié)點(diǎn)所對(duì)應(yīng)的殘值為

(5)

式中，R(e)為單元節(jié)點(diǎn)殘差，其他各參數(shù)意義與式(4)相同，利用形函數(shù)N對(duì)單元內(nèi)的場(chǎng)變量e進(jìn)行插值

e=N·en

(6)

式中，en為單元節(jié)點(diǎn)處的能量密度，利用微分運(yùn)算法則對(duì)式(5)等號(hào)右側(cè)第一項(xiàng)進(jìn)行變換

(7)

由散度定理，式(7)等號(hào)右側(cè)第一項(xiàng)可表為

(8)

從而式(5)可表為

∫∫ΩηωNTNdA·en-∫∫ΩNTπin(x,y)dA

(9)

令應(yīng)變矩陣為

BT=

(10)

則表征薄板單元能量傳遞效應(yīng)的單元矩陣為

(11)

表征薄板單元能量耗散效應(yīng)的單元矩陣為

(12)

令

(13)

與外界能量輸入有關(guān)的單元輸入功率向量為

F(e)=∫∫ΩNTπin(x,y)dA

(14)

與單元邊界能量密度梯度有關(guān)的向量為

(15)

從而薄板單元的節(jié)點(diǎn)殘差可表為

R(e)=K(e)en-Q(e)-F(e)

(16)

在計(jì)算得到各單元的有限元方程后，可按照各單元對(duì)節(jié)點(diǎn)殘值的貢獻(xiàn)關(guān)系將各單元所對(duì)應(yīng)的式(16)進(jìn)行組裝，其過(guò)程與靜力學(xué)有限元求解的有限元方程組裝相似，在組裝時(shí)由于能量密度為不具有方向性的一維場(chǎng)變量，且在計(jì)算單元矩陣時(shí)取各個(gè)單元的局部坐標(biāo)系坐標(biāo)軸相互平行，因此不需再進(jìn)行坐標(biāo)變換。與單元邊界能量密度梯度有關(guān)的向量Q(e)反映了經(jīng)由單元邊界傳導(dǎo)的功率流，容易看出，在進(jìn)行單元組裝，建立薄板的整體有限元方程式的過(guò)程中，位于薄板內(nèi)部單元的Q(e)將相互抵消而不對(duì)薄板的能量密度分布產(chǎn)生影響，組裝得到的整體有限元方程為

R=Ke-Q-F

(17)

式中：K為對(duì)各單元矩陣進(jìn)行組裝后得到的整體矩陣；F為外界對(duì)薄板的功率輸入向量；e為節(jié)點(diǎn)能量密度向量；若薄板的邊界無(wú)能量流出，則邊界處的Q(e)亦為0。為使得所建立的整體有限元方程盡可能滿足薄板能量密度控制方程，各節(jié)點(diǎn)殘值應(yīng)為0，從而薄板的整體有限元方程為

Ke=F

(18)

2 任意形狀板的能量密度求解

2.1 線性四邊形單元

功率流有限元方法視能量密度基本未知量，不再關(guān)注于結(jié)構(gòu)振動(dòng)的模態(tài)與波長(zhǎng)，因此網(wǎng)格劃分的密度不再受激勵(lì)頻率的影響，當(dāng)輸入的頻率增大時(shí)不需再劃分過(guò)于細(xì)密的網(wǎng)格。薄板能量密度為無(wú)方向性的一維場(chǎng)變量，在使用四節(jié)點(diǎn)的線性四邊形單元時(shí)，其單元形函數(shù)的形式為

N=[N1N2N3N4]

(19)

如圖3所示，可將劃分得到的物理坐標(biāo)下的四邊形單元映射為自然坐標(biāo)下的等參單元。

圖3 薄板等參單元的轉(zhuǎn)換Fig. 3 Thin plate’s isoparametric element transition

在自然坐標(biāo)下線性四邊形單元的形函數(shù)可表為自然坐標(biāo)(ξ,ζ)的函數(shù)，若第m個(gè)(m=1, 2, 3, 4)節(jié)點(diǎn)的自然坐標(biāo)為(ξm,ζm)，則該點(diǎn)所對(duì)應(yīng)的形函數(shù)為

(20)

物理坐標(biāo)系與自然坐標(biāo)系下面積微元之間的關(guān)系為

(21)

J為坐標(biāo)映射雅可比矩陣

J=[?ξN?ζN]TX

(22)

式中，X為線性四邊形單元節(jié)點(diǎn)的物理坐標(biāo)以行排列的形式組成的矩陣，則應(yīng)變矩陣

(23)

在將四邊形單元由物理坐標(biāo)映射到自然坐標(biāo)后，相應(yīng)的單元矩陣可在自然坐標(biāo)系下計(jì)算得到，其中表征單元能量傳遞特性的矩陣為

(24)

單元的能量耗散矩陣為

(25)

式(24)與式(25)中的積分均可借助于數(shù)值解法求解。當(dāng)薄板為規(guī)則的矩形時(shí)，使用規(guī)則的矩形單元便可實(shí)現(xiàn)網(wǎng)格劃分，此時(shí)式(24)與式(25)中的積分為多項(xiàng)式積分，借助高斯積分可計(jì)算得到精確的單元矩陣，通過(guò)單元矩陣的組裝求解便可得到所劃分的矩形單元節(jié)點(diǎn)處能量密度響應(yīng)。四邊形單元亦可應(yīng)用于其他不規(guī)則求解域，但此時(shí)便難以再劃分出矩形單元。當(dāng)采用非矩形的四邊形單元時(shí)，單元所對(duì)應(yīng)的雅可比矩陣不再為常數(shù)陣，式(24)中的積分也不再為多項(xiàng)式積分，因此使用高斯積分不能得到式(24)精確值。但當(dāng)所劃分出的四邊形單元與矩形單元較為接近時(shí)，可以認(rèn)為雅可比矩陣在單元內(nèi)的變化較小，在將其視為常數(shù)矩陣的情況下仍可使用高斯積分進(jìn)行計(jì)算。

以一不規(guī)則四邊形板為計(jì)算實(shí)例，該四邊形板的左右兩邊長(zhǎng)度分別為1 m，1.12 m，上下兩邊長(zhǎng)度分別為1.5 m，1 m，由于幾何形狀較為規(guī)則，使用四邊形單元對(duì)其進(jìn)行網(wǎng)格劃分。計(jì)算所使用參數(shù)如表1。

表1 薄板參數(shù)表

在該四邊形板的形心位置處模擬外界激勵(lì)施加集中單位功率激勵(lì)，并使用δ函數(shù)表示該形式的激勵(lì)

πinP·δ(x-x0)·δ(y-y0)

(26)

從而對(duì)于存在外界功率激勵(lì)的單元，與外界能量輸入有關(guān)的單元輸入功率向量為

F(e)=∫∫ΩN(ξ0,ζ0)T·P·δ(ξ-ξ0)·δ(ζ-ζ0)dA=

P·[N(ξ0，ζ0)]T

(27)

在式(24)中，薄板的彎曲波群速度cg為

(28)

式中，D為薄板的彎曲剛度。

使用高斯積分別對(duì)式(24)和式(25)進(jìn)行求解，在得到各個(gè)單元的有限元矩陣后通過(guò)組裝與求解可得到節(jié)點(diǎn)處的能量密度響應(yīng)，進(jìn)一步的利用式(6)進(jìn)行插值可得單元內(nèi)部任意位置處的能量密度。為更直觀的對(duì)比各位置處能量密度響應(yīng)大小，可采用功率級(jí)表示能量密度的大小

LP=10·lg(P/P0)

(29)

式中，P0為基準(zhǔn)能量密度，取P0=1×10-12J/m2。在平面波遠(yuǎn)場(chǎng)疊加假設(shè)下，經(jīng)過(guò)周期平均以及局部空間平均后的能量密度e與功率流q之間關(guān)系為

(30)

也即功率流沿能量密度的負(fù)梯度方向，將式(6)代入式(30)，可得

(31)

圖4為薄板在單點(diǎn)激勵(lì)下的能量密度響應(yīng)以及功率流分布情況，由圖4可知，能量密度響應(yīng)在激勵(lì)點(diǎn)處取得最大值，由于結(jié)構(gòu)阻尼的效應(yīng)，能量經(jīng)由功率流傳遞的過(guò)程中不斷損耗，能量密度隨著與激勵(lì)點(diǎn)距離的增加而不斷衰減。功率流以激勵(lì)位置為源頭，流向結(jié)構(gòu)內(nèi)部能量密度較低的區(qū)域。結(jié)構(gòu)振動(dòng)的這種能量傳遞機(jī)制類(lèi)似于熱傳遞形式，但同時(shí)也與熱傳遞形式存在一定區(qū)別，在熱傳導(dǎo)問(wèn)題中，熱能經(jīng)由熱對(duì)流等方式離開(kāi)結(jié)構(gòu)，結(jié)構(gòu)本身并不耗散熱能，但對(duì)于能量密度在板類(lèi)結(jié)構(gòu)中的分布問(wèn)題，結(jié)構(gòu)內(nèi)的能量除可經(jīng)過(guò)聲輻射以及與其他結(jié)構(gòu)的耦合向外傳遞外，結(jié)構(gòu)阻尼亦耗散部分能量，當(dāng)不考慮單板結(jié)構(gòu)的能量流出時(shí)，外界的能量輸入與阻尼耗散引起的能量損耗之間形成平衡。

圖4 異形求解域的能量密度分布與能流Fig. 4 Energy density and flux of the girregular domain

2.2 線性三角形網(wǎng)格

更為復(fù)雜的求解域難以獲得質(zhì)量較高的四邊形網(wǎng)格，若此時(shí)仍采用四邊形網(wǎng)格，單元畸形的發(fā)生將使得計(jì)算結(jié)果誤差較大，此時(shí)應(yīng)該引入適應(yīng)性更強(qiáng)的三角形網(wǎng)格。

對(duì)于三節(jié)點(diǎn)的線性三角形單元，其形函數(shù)的形式為

N=[N1N2N3]

(32)

此時(shí)可取各節(jié)點(diǎn)的形函數(shù)為其對(duì)應(yīng)的面積坐標(biāo)，如圖5所示，也即

(33)

式中，Ae為三角形單元面積。

(34)

圖5 薄板三角形單元與面積坐標(biāo)Fig. 5 Triangle element and area coordinate

各節(jié)點(diǎn)所對(duì)應(yīng)的面積亦可以類(lèi)似求得。物理坐標(biāo)下線性三角形單元的能量耗散矩陣為

(35)

利用面積坐標(biāo)在三角形中的積分公式

(36)

式中，Lk為三角形單元的第k個(gè)節(jié)點(diǎn)所對(duì)應(yīng)的面積坐標(biāo)，可得

(37)

線性三角形單元的能量擴(kuò)散矩陣為

(38)

線性三角形單元的應(yīng)變矩陣為常數(shù)矩陣，從而式(38)為常數(shù)積分，將式(33)代入式(11)并在三角形單元內(nèi)進(jìn)行積分可得

(39)

其中，

(40)

(41)

(42)

式中，i,j,k取值為1, 2, 3，輪流換位，正序排列。

為驗(yàn)證三角形網(wǎng)格的適用性，再次使用三角形網(wǎng)格對(duì)圖4中的薄板進(jìn)行網(wǎng)格劃分求解并與四邊形網(wǎng)格以及統(tǒng)計(jì)能量法的求解結(jié)果進(jìn)行對(duì)比。仍采用表1中的薄板參數(shù)，計(jì)算頻率為1～5 kHz的1/3倍頻程中心頻率。由于統(tǒng)計(jì)能量法基于模態(tài)原理，提供的是子系統(tǒng)的平均能量響應(yīng)，而功率流有限元方法的建立基于波動(dòng)原理，得到的是具有空間分布特性的能量密度，為實(shí)現(xiàn)二者的對(duì)比，在將板結(jié)構(gòu)視為單一子系統(tǒng)的情況下可將具有空間分布特性的能量密度在子系統(tǒng)內(nèi)進(jìn)行空間平均

(43)

將式(6)代入式(43)可得

(44)

式中：M為網(wǎng)格數(shù)目；Ni和ei分別為第i個(gè)單元的形函數(shù)以及單元節(jié)點(diǎn)能量密度值。圖6為三者的對(duì)比結(jié)果，由圖可見(jiàn)，使用三角形單元以及四邊形單元時(shí)，功率流有限元空間平均后的能量密度皆與SEA結(jié)果吻合良好，表明三角形單元以及四邊形單元在功率流有限元中具有較高的精度以及良好的收斂性。

圖6 功率流有限元與統(tǒng)計(jì)能量法結(jié)果對(duì)比Fig.6 Comparison between results of PFFEA and SEA

三角形單元對(duì)復(fù)雜求解域具有更好的適應(yīng)性，可對(duì)較為復(fù)雜的求解域進(jìn)行網(wǎng)格離散，因此當(dāng)板類(lèi)結(jié)構(gòu)較為復(fù)雜時(shí)，可考慮采用三角形單元進(jìn)行網(wǎng)格劃分與求解。如圖7所示，以一類(lèi)似于我國(guó)內(nèi)陸輪廓形狀的薄板為實(shí)例，該薄板與前文所分析的方形薄板具有相同的材料屬性，由于其幾何形狀較為復(fù)雜，不適合使用四邊形單元對(duì)于求解域進(jìn)行離散，故選用三角形單元進(jìn)行網(wǎng)格劃分。

圖7 三角形網(wǎng)格的劃分(單位：m)Fig. 7 Triangular element mesh(unit：m)

利用式(37)和式(39)分別計(jì)算各個(gè)單元所對(duì)應(yīng)的能量耗散矩陣與能量擴(kuò)散矩陣，考慮無(wú)能量從該薄板流出的情形，則單元邊界能量密度梯度有關(guān)的向量Q(e)不需計(jì)算?？紤]多點(diǎn)激勵(lì)的形式，在該薄板上選取6處位置，施加單位集中功率激勵(lì)，與四邊形單元求解類(lèi)似，可利用δ函數(shù)表示集中功率激勵(lì)，進(jìn)而可求解得出外界對(duì)薄板功率輸入向量。通過(guò)對(duì)單元矩陣的組裝求解可得到節(jié)點(diǎn)處的能量密度響應(yīng)。

薄板內(nèi)能量密度分布如圖8所示，由圖可知，在集中功率輸入的位置能量密度較大，各個(gè)功率輸入點(diǎn)之間形成“能量密度峰”，在各個(gè)輸入點(diǎn)之間，隨著與功率激勵(lì)點(diǎn)位置距離的增加，能量在傳播過(guò)程中不斷耗散，能量密度也隨之衰減，從而在各個(gè)“能量密度峰”之間形成“能量密度谷”。

圖9為該薄板能量密度分布的等高線以及功率流分布圖，在多點(diǎn)激勵(lì)情況下，板內(nèi)的能量密度以激勵(lì)位置為圓心呈輻射狀分布。圖9中左側(cè)的兩處激勵(lì)位置雖與其余四處激勵(lì)位置具有相同的激勵(lì)，但由于右側(cè)兩處功率激勵(lì)位置處能量傳遞影響范圍較大，因此該激勵(lì)點(diǎn)位置能量分布較為稀疏，激勵(lì)點(diǎn)周?chē)芰棵芏容^小，右側(cè)四處激勵(lì)點(diǎn)由于能量傳遞影響范圍較小且激勵(lì)點(diǎn)相對(duì)集中，能量密度分布較為密集，因此激勵(lì)點(diǎn)附近能量密度較大。在多點(diǎn)激勵(lì)下，各個(gè)激勵(lì)點(diǎn)為功率流的源頭，功率流由各個(gè)激勵(lì)點(diǎn)流向板內(nèi)能量密度較低的區(qū)域。

薄板內(nèi)的能量密度由形函數(shù)直接插值得到，而不再類(lèi)似于靜力學(xué)有限元求解中需使用應(yīng)變矩陣以及本構(gòu)關(guān)系得到單元內(nèi)的應(yīng)力分布，三角形單元形函數(shù)的連續(xù)性保證了單元內(nèi)能量密度分布的連續(xù)性，因此不同于三角形單元在靜力學(xué)有限元求解時(shí)的“常應(yīng)力單元”現(xiàn)象，功率流有限元中所引入的三角形單元不是常能量密度單元。同時(shí)，形函數(shù)的德?tīng)査瘮?shù)性質(zhì)和單位分解性質(zhì)使得能量密度在不同單元的邊界保持連續(xù)。當(dāng)使用三角形單元時(shí)，同樣可由式(31)可確定板內(nèi)的功率流分布情況，此時(shí)值得注意的是，線性三角形單元的形函數(shù)為關(guān)于坐標(biāo)的線性函數(shù)，其單元應(yīng)變矩陣為常數(shù)矩陣，因此在三角形單元雖不為常能量密度單元，但卻為常功率流單元，當(dāng)需要得到較為精確的功率流分布時(shí)，應(yīng)增加三角形單元數(shù)目，或使用高次三角形單元以及四邊形單元等單元類(lèi)型。

3 結(jié) 論

利用功率流有限元方法對(duì)結(jié)構(gòu)的振動(dòng)響應(yīng)進(jìn)行分析具有諸多優(yōu)勢(shì)，本文在已有研究的基礎(chǔ)上，利用該方法對(duì)薄板在功率激勵(lì)下的能量密度的分布進(jìn)行了研究，分別使用四邊形單元和三角形單元對(duì)薄板結(jié)構(gòu)進(jìn)行網(wǎng)格劃分，建立了對(duì)應(yīng)的單元矩陣和有限元方程，求解有限元方程得出了薄板的能量密度響應(yīng)。

圖8 薄板內(nèi)的能量分布Fig. 8 Energy density distribution of the thin plate

圖9 薄板的能量密度分布與能流Fig. 9 Energy density and energy flux of the thin plate

四邊形單元可對(duì)應(yīng)用廣泛的矩形薄板進(jìn)行網(wǎng)格劃分求解，同時(shí)當(dāng)求解域非矩形但劃分得到的四邊形單元畸形程度較低時(shí)仍可使用四邊形單元。當(dāng)求解域較為復(fù)雜時(shí)，應(yīng)使用三角形單元對(duì)薄板進(jìn)行網(wǎng)格劃分時(shí)。三角形單元的引入可將功率流有限元法拓展至任意形狀的薄板類(lèi)結(jié)構(gòu)，同時(shí)三角形單元內(nèi)的能量密度由形函數(shù)插值得到，因此在功率流有限元中利用三角形單元得到的能量密度連續(xù)分布，并不出現(xiàn)“常能量密度單元”。三角形單元為常功率流單元，因此當(dāng)需要得到較為詳細(xì)的功率流分布時(shí)，應(yīng)使用較為細(xì)密的三角形網(wǎng)格或采用其他單元類(lèi)型。多種單元類(lèi)型的使用可有效的將功率流有限元方法應(yīng)用于異形薄板以及復(fù)雜耦合板結(jié)構(gòu)。

[ 1 ] 宋孔杰，張蔚波，牛軍川. 功率流理論在柔性振動(dòng)控制技術(shù)中的應(yīng)用與發(fā)展 [J]. 機(jī)械工程學(xué)報(bào), 2003, 39(9): 23-28. SONG Kongjie, ZHANG Weibo, NIU Junchuan. Application and development of power flow theories in the field of the vibration control for flexible system [J]. Chinese Journal of Mechanical Engineering, 2003, 39(9): 23-28.

[ 2 ] CHO P E. Energy flow analysis of coupled structures [D]. United States: Purdue University, 1993.

[ 3 ] 朱繼清. 基于能量流有限元的耦合結(jié)構(gòu)振動(dòng)傳遞特性分析 [D]. 濟(jì)南：山東大學(xué), 2011.

[ 4 ] NEFSKE D J, SUNG S H. Power flow finite element analysis of dynamic systems: basic theory and application to beams [J]. Journal of Vibration & Acoustics, 1989, 111(1): 94-100.

[ 5 ] WOHLEVER J C, BERNHARD R J. Mechanical energy flow models of rods and beams [J]. Journal of Sound and Vibration, 1992, 153(1): 1-19.

[ 6 ] CHO P E, BERNHARD R J. Energy flow analysis of coupled beams [J]. Journal of Sound and Vibration, 1998, 211(4): 593-605.

[ 7 ] SONG J H, HONG S Y, KANG Y, et al. Vibrational energy flow analysis of penetration beam-plate coupled structures [J]. Journal of Mechanical Science and Technology, 2011, 25(3): 567-576.

[ 8 ] BOUTHIER O M, BERNHARD R J. Simple models of energy flow in vibrating membranes [J]. Journal of Sound and Vibration. 1995, 182(1): 129-147.

[ 9 ] BOUTHIER O M, BERNHARD R J. Simple models of the energetics of transversely vibrating plates [J]. Journal of Sound and Vibration, 1995, 182(1): 149-166.

[10] PARK D H, HONG S Y, KIL H G, et al. Power flow models and analysis of in-plane waves in finite coupled thin plates [J]. Journal of Sound and Vibration, 2001, 244(4): 651-668.

[11] 李坤朋. 基于能量有限元法的板耦合結(jié)構(gòu)振動(dòng)特性分析 [D]. 濟(jì)南：山東大學(xué), 2013.

[12] NIU Junchuan, LI Kunpeng. Energy flow finite element analysis of L-shaped plate including three types of waves [J]. Applied Mechanics and Materials, 2013, 353/356: 3365-3368.

[13] NIU Junchuan, LI Kunpeng. Energy finite element analysis of n-shaped plate structures with three types of wave [J], Journal of Vibration Engineering and Technologies, 2015, 3(5): 615-625.

[14] 江民圣，牛軍川，鄭建華，等. L 型耦合板結(jié)構(gòu)能量傳遞系數(shù)特性的研究 [J]. 振動(dòng)與沖擊, 2015,34(17): 131-136.JIANG Minsheng, NIU Junchuan, ZHENG Jianhua, et al. Energy transfer coefficients features of L-shape coupled plates [J]. Journal of Vibration and Shocks, 2015, 34(17): 131-136.

[15] SEO S, HONG S, KIL H. Power flow analysis of reinforced beam-plate coupled structures [J]. Journal of Sound and Vibration. 2003, 259(5): 1109-1129.

[16] PEREIRA V S, DOS SANTOS J M C. Coupled plate energy models at mid-and high-frequency vibrations [J]. Computers & Structures, 2014, 134(4): 48-61.

Energy density analysis of irregular shaped plates based on the power flow finite element method

LIU Zhihui1, NIU Junchuan1,2, ZHOU Yiqun1

(1. School of Mechanical Engineering, Shandong University, Ji’nan 250061, China; 2. Key Laboratory of High Efficiency and Clean Mechanical Manufacture, Shandong University, Ji’nan 250061, China)

It has important significance to analyze and predict the vibration response of structures under different frequencies. The power flow finite element method has become a research focus for vibration analysis due to the advantages of broad applicable frequency range and the detailed information provided. The power flow finite element method was implemented to solve the flexural wave energy density of a thin plate, and the weighted residual method was used to derive the residual of the node points. The linear quadrilateral mesh was used to partition the thin plate and the element’s finite element equation was derived, then the global finite element equation was assembled and solved, and the energy density of the nodes was obtained. The linear triangular element was introduced to partition the plate with complex shape. The research is meaningful for the implementation of the power flow element method on thin plates with arbitrary shapes.

vibration; power flow finite element analysis; irregular shaped plate; energy density

國(guó)家自然科學(xué)基金(51275275；51675306)；山東省優(yōu)秀中青年科學(xué)家科研獎(jiǎng)勵(lì)基金(BS2010ZZ006)

2016-01-19 修改稿收到日期： 2016-06-17

劉知輝 男，碩士，1992年7月生

牛軍川 男，教授，博士生導(dǎo)師，1974年生

TH212；TH213.3

A

10.13465/j.cnki.jvs.2017.16.029