于開平 王懷志

(哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 150001)

0 引言

隨著技術(shù)的發(fā)展，高速行駛的機車以及航空航天中高速飛行器等結(jié)構(gòu)的研究都對動響應(yīng)預(yù)示方法提出了更多的需求，傳統(tǒng)振動環(huán)境預(yù)示方法，如有限元和邊界元等[1]，隨著頻率的提高需要更高的網(wǎng)格密度，導(dǎo)致單元類的方法在中高頻響應(yīng)計算上效率較低，且隨著頻率的提高，結(jié)構(gòu)的模態(tài)更加密集，有限元方法很難得到準(zhǔn)確的結(jié)構(gòu)固有頻率信息[2]。統(tǒng)計能量分析方法（SEA）的出現(xiàn)可以解決有限元等方法在中高頻響應(yīng)計算上的一些問題[3]，基于能量在頻率和空間上的平均，SEA 可以將系統(tǒng)簡化為有限個數(shù)的子系統(tǒng)，簡化后的模型自由度數(shù)與子系統(tǒng)數(shù)相關(guān)，相對于有限元方法，統(tǒng)計能量模型的自由度小的多，同時由于進行了空間和頻率上的取平均，SEA 也適用于具有一些不確定性和隨機特性的結(jié)構(gòu)[4]。但是由于子系統(tǒng)需要滿足模態(tài)密度的要求，子系統(tǒng)的特征尺寸一般較大，由于對結(jié)構(gòu)的能量進行空間平均，SEA 無法對得到結(jié)構(gòu)的局部響應(yīng)信息，也無法考慮結(jié)構(gòu)局部質(zhì)量、剛度、阻尼和載荷等的不均勻分布問題。

能量有限元方法（EFEA）最初由Nefske[5]通過對功率流方程的研究，得到了與熱傳導(dǎo)方程類似的能量密度控制方程，該方程可以使用有限元程序進行求解。Wohlever[6]推導(dǎo)了桿和梁結(jié)構(gòu)的能量密度方程，并對其特性進行了研究。Bouthier[7]推導(dǎo)了板的能量密度方程，并與有限元方法進行了比較。Bouthier[8]對板的能量有限元進行了系統(tǒng)的研究，并給出了無限大結(jié)構(gòu)的能量密度控制方程。Bouthier[9]推導(dǎo)了膜結(jié)構(gòu)的能量密度控制方程。Park[10]研究了基于Timoshenko 梁理論的 EFEA，HAN[11]將 Rayleigh-Love 和Rayleigh-Bishop 桿的結(jié)果與經(jīng)典桿理論的結(jié)果進行了比較，Bolt[12]用EFEA 對曲梁結(jié)構(gòu)進行了研究。Wang[13]考慮了熱作用下板的 EFEA，Chen[14]對梁在受到預(yù)應(yīng)力作用下的EFEA進行了研究，并給出了預(yù)應(yīng)力作用下修正的輸入功率就算公式。

EFEA 還被應(yīng)用于一些復(fù)雜的結(jié)構(gòu)，代文強等[15]對高速列車的車體和聲腔進行了建模，綜合考慮了車廂受到的輪軌激勵、二系懸掛力、輪軌噪聲以及氣動噪聲等復(fù)雜激勵的作用，與實驗結(jié)果進行了相比，在分析頻段內(nèi)仿真結(jié)果可以較好的預(yù)測車內(nèi)的聲壓級響應(yīng)，在文獻[16]中，他們進一步分析了各種激勵源對結(jié)構(gòu)和內(nèi)部聲場的影響，據(jù)此對輪軌噪聲進行了優(yōu)化，有效的降低了車體內(nèi)部的聲壓級。張雙雙[17]使用EFEA 對車內(nèi)噪聲進行了仿真分析。

在航空航天領(lǐng)域，Lima[18]使用EFEA 計算了飛機由湍流邊界層激勵引起的振動和內(nèi)部噪聲，建立了可用于EFEA 的湍流邊界層載荷模型，通過對實測飛行數(shù)據(jù)的比較驗證了EFEA 仿真結(jié)果的可靠性。Vlahopoulos[19]等使用EFEA 對鋁制圓柱殼進行了仿真分析，該結(jié)構(gòu)具有沿周向的周期性加筋和沿軸向的加筋，通過激振器在結(jié)構(gòu)上施加4 個單點激勵，與實際測試結(jié)果比較表明，EFEA 可以預(yù)測由加筋結(jié)構(gòu)帶來的周期特性，仿真結(jié)果和實測結(jié)果表現(xiàn)了良好的一致性，文中還比較了是否考慮加筋的周期效應(yīng)對功率傳輸系數(shù)的影響，結(jié)果表明，如果不考慮加筋的周期效應(yīng)，則傳輸?shù)哪芰繒兊酶?。Vlahopoulos[20]將EFEA 應(yīng)用于復(fù)合材料飛機機身結(jié)構(gòu)，其中復(fù)合材料使用了等效材料建模的方式，基于耦合邊界的位移協(xié)調(diào)和力平衡條件推導(dǎo)了復(fù)合材料層合板的耦合損耗傳輸系數(shù)，文中將飛機機身簡化為圓柱殼結(jié)構(gòu)，在10 個位置通過激振器施加單點隨機激勵，通過比較發(fā)現(xiàn)，EFEA 的仿真結(jié)果與實驗結(jié)果相關(guān)性很好，誤差基本滿足要求。

國內(nèi)對EFEA 在航空航天方面也有較為深入的應(yīng)用，Xie[21]使用殼單元建立了截錐形儀器艙的EFEA 模型，并分別對儀器艙在點載荷和風(fēng)洞實驗提供的壓力載荷作用下的聲振耦合特性進行了研究。王懷志、于開平等[22]首先針對雙星整流罩進行了EFEA 建模（圖1），該整流罩結(jié)構(gòu)具有兩個聲腔以及多個復(fù)合材料殼體結(jié)構(gòu)，實驗?zāi)M了整流罩處于高聲壓級作用的情況，整流罩外表面受到混響場聲激勵作用，通過與實驗值以及SEA 結(jié)果的比較，EFEA 在中高頻段獲得了良好的預(yù)示結(jié)果（圖2）。然后，他們又使用EFEA 對氣瓶整艙進行了詳細(xì)的EFEA 建模[23]，考慮了結(jié)構(gòu)受到外界混響聲場激勵作用下的聲振耦合特性，仿真結(jié)果與SEA 的結(jié)果以及實測結(jié)果進行了對比，驗證了EFEA 模型的可靠性。

圖1 整流罩結(jié)構(gòu)EFEA 模型及仿真結(jié)果Fig.1 EFEA model and simulation results of the fairing

圖2 聲腔的聲壓級響應(yīng)Fig.2 Sound pressure level response of the acoustic cavity

此外，原凱[24]、Xie[25]對EFEA 的發(fā)展以及應(yīng)用進行了總結(jié)。相對于SEA，EFEA 可以考慮更多結(jié)構(gòu)的局部特征（如幾何、材料屬性或者阻尼等的變化），并得到結(jié)構(gòu)的局部響應(yīng)。相對有限元方法，EFEA 不僅具有性能上的優(yōu)勢，而且對于中高頻分析，由于響應(yīng)對結(jié)構(gòu)參數(shù)變化比較敏感，有限元基于確定參數(shù)的建模往往會帶來較大的誤差[26-28]，基于統(tǒng)計的能量類方法可以較好的解決這個問題[3]。盡管EFEA 已經(jīng)有了較為深入的研究，但是對于EFEA 的一些基本假設(shè)，以及計算結(jié)果的可靠性缺乏系統(tǒng)的總結(jié)，本文以板為例，分析了目前在EFEA 中使用的一些基本假設(shè)，并總結(jié)了目前文獻中經(jīng)過驗證的一些可靠性評估方法，對于不滿足平面波假設(shè)的情況，對一種改進方法進行了參數(shù)研究以及仿真分析。

1 能量密度的控制方程

EFEA 的能量密度控制方程基于功率流平衡，這種平衡關(guān)系可以用單元體來表示，如圖3 所示。

圖3 單元體的能量流關(guān)系Fig.3 Energy flow relationship of unit element

這個關(guān)系可以用以下公式來表示[5]

其中E表示單元體的能量，表示單元體邊界流出的能量強度，πdiss表示單元體內(nèi)由于阻尼等耗散的能量，πin表示單元輸出的能量。

假設(shè)系統(tǒng)阻尼可以簡化為滯后阻尼，且阻尼足夠小，則耗散的能量可以用以下公式來表述[3]

其中，η是內(nèi)損耗因子，ω是圓頻率。

結(jié)構(gòu)中單個面波分量可以用以下公式來描述

其中，A為波的幅值，K表示波的波數(shù)，r表示距離原點的距離。

以彎曲波為例，能量密度和能量強度可以表述為

其中，D表示彎曲剛度，ρ表示密度，k1為波數(shù)K的實部。

當(dāng)阻尼很小時，波數(shù)K的虛部可以忽略，則通過公式(4)可以得到能量密度和能量強度之間的關(guān)系為

其中，Cg是彎曲波對應(yīng)的群速度。

將公式(5)代入公式(1)中，可以得到單彎曲面波的能量密度控制方程

僅考慮穩(wěn)態(tài)振動的情況，可以去除上式中的時間微分項，則穩(wěn)態(tài)的能量密度控制方程可以表述為

可以發(fā)現(xiàn)，從公式(1)到(7)使用面波方法的推導(dǎo)并不涉及具體的結(jié)構(gòu)（如桿、梁和板等），這說明公式(6)和(7)適用于普遍的擴散場面波的傳播過程。面波假設(shè)廣泛應(yīng)用于能量類方法的公式推導(dǎo)中[4]，EFEA 假設(shè)結(jié)構(gòu)內(nèi)混響場占主導(dǎo)，可以近似用面波來表示。

如果將考慮反射波[6]或者將結(jié)構(gòu)擴展為多維結(jié)構(gòu)[7]，在能量密度和能量強度的表達(dá)式中將會出現(xiàn)耦合項，該耦合項可以通過進行局部空間平均的方法來消除掉。所以，理論上，EFEA 與SEA類似，是需要進行空間平均的，不同的是EFEA的空間平均是局部平均，最小需要在一個波長內(nèi)進行平均即可。對于桿、歐拉-伯努力梁[6]、板[7]、圓柱殼[29]、鐵木辛柯梁[10]以及一階剪切變形板[30]（單波形式）等常見結(jié)構(gòu)，通用的能量密度控制方程可以表述為

其中，〈*〉表示需要對該值進行一個波長上的空間平均。而公式(5)表示的能量強度則可以重新寫成為

能量密度方程可以用有限元方法進行離散，離散方法可以使用常規(guī)的EFEA，也可以使用0階EFEA[31]。0 階EFEA 采用了類似數(shù)值離散的集中質(zhì)量方法，每個單元可以看作是一個小的SEA 子系統(tǒng)，該方法可以直接使用耦合損耗因子，在相同網(wǎng)格密度的情況下，等效剛度矩陣由于是對角陣，實際上計算量略小于傳統(tǒng)的有限元離散方法。

公式(8)可通過Galerkin 加權(quán)余量方法來構(gòu)建能量有限元方程，對于2 維板結(jié)構(gòu)，余量方程為

其中N i為型函數(shù)。單元的能量可以表示為

公式(10)中的第一個積分項可以展開為

其中，

通過展開可以將公式(10)可以轉(zhuǎn)化為

公式(14)可以寫成以下形式

其中，

Q e表示單元邊界上輸入的功率流，對于結(jié)構(gòu)內(nèi)部（不存在幾何、材料等不連續(xù)），該項可以直接消去。而對于耦合邊界，需要通過耦合邊界的功率流平衡關(guān)系來得到耦合傳遞矩陣[32]，則板的能量有限元公式可以表示為

其中，JC表示耦合邊界的能量傳輸矩陣。

對于公式中的能量輸入項Pin，一般可以通過無限大結(jié)構(gòu)的原點輸入導(dǎo)納來得到，即

如果使用相同無限大結(jié)構(gòu)的原點輸入導(dǎo)納計算輸入功率，EFEA 的計算結(jié)果實際相當(dāng)于SEA進行頻率平均后的結(jié)果[3]。因為對于有限大結(jié)構(gòu)，如果在一定頻帶內(nèi)的，載荷為隨機或者系統(tǒng)本身為隨機，則其在帶寬內(nèi)的平均輸入導(dǎo)納約等于無限大板的原點輸入導(dǎo)納。這是EFEA 需要進行頻率平均的一個原因。

2 EFEA 結(jié)果的有效性

目前，EFEA 計算結(jié)果的有效性主要可以從以下幾個方面：能量耗散的計算誤差，特征長度與結(jié)構(gòu)內(nèi)最大作用波長的比值以及直接場能量的占比。

Wohlever[6]指出，在EFEA 中，能量密度是結(jié)構(gòu)動能和勢能的和，而損耗的能量主要來自結(jié)構(gòu)的振動，因此，公式中的能量密度實際代表的是結(jié)構(gòu)動能密度的2 倍，對于結(jié)構(gòu)動能和勢能同相位的結(jié)構(gòu)，如梁和板，在小阻尼假設(shè)下，可以近似認(rèn)為動能和勢能是處處相等的，但是對于桿結(jié)構(gòu)，穩(wěn)態(tài)振動時動能和勢能具有相反的相位，所以如果桿結(jié)構(gòu)內(nèi)含有波長數(shù)不是半波長的整數(shù)倍，能量耗散值的計算就會引入誤差。當(dāng)頻率比較高時，桿結(jié)構(gòu)中含有的波長數(shù)增加，非整數(shù)波長引入的誤差比例會逐漸降低。

Gur[33]通過對板梁結(jié)構(gòu)在特定頻帶內(nèi)響應(yīng)的研究提出了板和梁結(jié)構(gòu)的波長準(zhǔn)則，對于梁結(jié)構(gòu)，一般要求單個梁的長度至少大于5 倍的彎曲波長；對于板結(jié)構(gòu)，板的特征長度要求大于2.4 倍的彎曲波長。

Moens[34]通過對不同形狀的板的研究指出，如果板結(jié)構(gòu)的特征長度小于2.4，動能和勢能的相等性就無法得到保證。

此外，Langley[35]指出，在強阻尼情況下，直接場在結(jié)構(gòu)響應(yīng)中占比很大，在板結(jié)構(gòu)中，由EFEA能量密度方程描述的能量密度分布是成正比，而實際的能量密度分布應(yīng)該是成正比，其中r柱坐標(biāo)中響應(yīng)點與激勵點之間的距離。因此，在直接場占比較大的情況下，如關(guān)注頻率較高時或者阻尼比較大時，相對于實際值EFEA 的結(jié)果一般表現(xiàn)為，在激勵點附近較小而在離激勵點很遠(yuǎn)的地方則相對較大。

Xiang jie[36]通過輻射能量傳輸方程給出了驗證板結(jié)構(gòu)EFEA 結(jié)果有效性的準(zhǔn)則（公式(19)），文中指出EFEA 基于混響場假設(shè)，而公式(19)可以用來判斷結(jié)構(gòu)是否滿足混響場假設(shè)。對于阻尼越大，頻率越高，板的特征尺寸越大，直接場部分耗散的能量會更多，而進入混響場的能量就會相對變小，從而導(dǎo)致不滿足文中提出的混響場假設(shè)

3 非平面波假設(shè)的改進

第2 節(jié)的推導(dǎo)方法是EFEA 常用方法，該方法實際需要假設(shè)結(jié)構(gòu)內(nèi)的波滿足平面波假設(shè)[7,8]。根據(jù)文獻[4]和[36]，在點激勵作用下，如果直接擴散場占主導(dǎo)，則平面波假設(shè)是不滿足的，而混響場占主導(dǎo)時，由于各個方向的波是均勻分布的，可以滿足平面波假設(shè)。因此，可以考慮將板的響應(yīng)分為直接擴散場和混響場兩個部分[37,38]，即

其中，下標(biāo)“total”表示總能量，“dir”代表直接場的能量，“rev”代表混響場的能量。

直接擴散場的能量密度可以用公式(21)來表述[8]，而混響場部分的能量可以通過公式(8)來描述。

其中，r為距離激勵點的距離。

以一個四邊簡支單板為例，板的尺寸為1.5m×1.5m，厚度為0.002 m，材料參數(shù)如下：密度2700kg/m3，楊氏模量71×109N/m2，泊松比為0.33。在中心點位置施加一個頻帶400-500 Hz 區(qū)間的隨機點激勵，力的幅值大小為10 N。在400-500 Hz 區(qū)間內(nèi)，板共有33 階模態(tài)，有限元能量結(jié)果在頻帶內(nèi)取平均值，然后將結(jié)果在一個波長上取平均值。

對于單板系統(tǒng)，直接場能量在到達(dá)結(jié)構(gòu)邊界的時候會發(fā)生反射，而反射的功率流則作為混響場的能量輸入，即在結(jié)構(gòu)的邊界，能量流滿足公式(22)

板在x=0.75 m 中線上的響應(yīng)如圖4 所示。圖注中的“averaged”表示有限元平均后的結(jié)果，“rev”表示EFEA 計算的混響場能量，“dir”表示公式(21)計算的直接場能量，實心○標(biāo)注的是基于面波公式的結(jié)果。圖4 結(jié)果所采用的滯后阻尼系數(shù)為0.01，根據(jù)公式(19)計算得到判定系數(shù)為0.32，值略小于1。可以發(fā)現(xiàn)，此時直接場的能量相對于混響場的能量小了很多，與有限元結(jié)果相比，EFEA 基于面波的公式計算結(jié)果在激勵點位置偏小，在邊界位置偏大，這與上節(jié)的預(yù)期相符[8,35]。同時由于在基于面波的公式推導(dǎo)中忽略了近場項，這也導(dǎo)致了EFEA 結(jié)果在激勵點和邊界位置出現(xiàn)一定的偏差。此外，有限大板在中線和對角線位置的響應(yīng)要比其他位置大，所以用中線位置的結(jié)果進行對比看起來差別較大，實際上差別較大的位置主要集中在激勵點附近以及中線和對角線位置，其他位置差別并不是很大，各種方法得到板的整體能量是基本相同的。整體結(jié)果對比可以得到，EFEA 基于面波的結(jié)果與SEA計算的能量密度平均值基本差別不大，而考慮了直接場能量的結(jié)果與有限元結(jié)果整體吻合較好。

圖4 板在x=0.5 m 中線位置的響應(yīng)（滯后阻尼系數(shù)為0.01）Fig.4 Response of the plate at x=0.5 m with damp loss factor is 0.01

將阻尼系數(shù)設(shè)置為0.1，此時由公式(19)計算的判定系數(shù)為3.19，值略大于1，結(jié)果如圖5 所示。可以發(fā)現(xiàn)，此時混響場的能量相對直接場的能量整體小了很多，兩者疊加后的總能量與有限元結(jié)果差別較小，基于面波的公式此時仍然表現(xiàn)為在激勵點處值偏小而在邊界位置值偏大，但與SEA 給出的平均值比較，EFEA 給出的能量密度分布相對更接近有限元的結(jié)果。

圖5 滯后阻尼系數(shù)為0.1 時，板在x=0.5 m 中線位置的響應(yīng)Fig.5 Response of the plate at x=0.5 m with damp loss factor 0.1

通過兩組結(jié)果的比較可以發(fā)現(xiàn)，采用直接場疊加混響場能量的結(jié)果整體與有限元的結(jié)果更為吻合，同時也可以發(fā)現(xiàn)，滯后阻尼系數(shù)為0.01時，誤差相對更大一些，這是由于在EFEA 公式的推導(dǎo)中，一般會忽略近場項，當(dāng)阻尼系數(shù)較大時，近場項很快就衰減到很小值，而阻尼系數(shù)較小時，近場項的影響就更加明顯。公式(21)給出的結(jié)果實際上僅當(dāng)k·r遠(yuǎn)大于1 時是有效的[39]，公式在r=0的位置是奇異的，因此在k·r接近于1 的值附近，可以使用r=0 附近的平均值來代替，在本例中，r=0 附近的值用k·r=2 的平均值來代替。

類似的，對圖6 所示的四邊簡支耦合板進行能量有限元分析，兩個板的尺寸均為1.5m×1.5m，厚度分別為2mm 和4mm，材料參數(shù)如下：密度2700kg/m3，楊氏模量71×109N/m2，泊松比為0.33。在其中一個板的中心位置(圖6 所示)施加一個作用于400-500Hz 區(qū)間內(nèi)隨機點激勵，通過有限元分析可得耦合模型在400-500Hz 區(qū)間內(nèi)共有49階模態(tài)。

圖6 四邊簡支耦合板Fig.6 Coupling plates with simply supported on all edges

耦合板在內(nèi)部邊界以及耦合邊界滿足以下能量流平衡條件，

假設(shè)直接場傳入其他板的能量滿足平面波假設(shè)，則當(dāng)滯后阻尼系數(shù)為0.01 時，耦合板的在x=0.5 m 位置的響應(yīng)如圖7 所示。公式(19)計算的判定系數(shù)和單板相同，此時直接場能量除了在激勵點位置附近外，其他位置的能量相對小很多，所以基于面波的公式在距離激勵點較遠(yuǎn)的位置（例如靠近邊界以及其他板）誤差相對較小。

圖7 滯后阻尼系數(shù)為0.01 時，耦合板在x=0.5 m 中線位置的響應(yīng)Fig.7 Response of the plate at x=0.5 m with damp loss factor 0.01

圖8 滯后阻尼系數(shù)為0.1 時，耦合板在x=0.5 m 中線位置的響應(yīng)Fig.8 Response of the plate at x=0.5 m with damp loss factor 0.1

滯后阻尼系數(shù)為0.1 時，直接場的能量占主導(dǎo)，此時在遠(yuǎn)離激勵點位置的能量以及混響場的能量幾乎可以忽略?；谄矫娌ǖ墓皆谶@種情況下的結(jié)果相對SEA 具有一定的優(yōu)勢。對于耦合板，不同阻尼的情況下，通過單獨計算直接場能量和混響場能量的方法結(jié)果與有限元結(jié)果更加接近，這與四邊耦合單板的結(jié)果時一致的。

4 結(jié)論

本文簡述了能量有限元研究進展，并通過平面波假設(shè)，說明了以波動解為基礎(chǔ)的物理模型都可以使用能量有限元的推導(dǎo)方法得到能量平衡關(guān)系。對于目前在EFEA 中使用的一些可靠性評估方法進行了系統(tǒng)的解釋。對于非平面波的情況，如點載荷作用下的板，以四邊簡支單板和耦合板為例，對板在直接場和混響場占主兩種導(dǎo)情況下（分別對應(yīng)兩種阻尼大?。┑捻憫?yīng)結(jié)果進行對比分析，計算結(jié)果與具有足夠網(wǎng)格密度的有限元結(jié)果進行了比較，結(jié)果表明，對于單點載荷作用下的平板結(jié)構(gòu)，直接場的能量對激勵點附近的響應(yīng)結(jié)果影響很大，而阻尼足夠大時，結(jié)構(gòu)中的能量以直接場的能量占主導(dǎo)，傳統(tǒng)平面波假設(shè)的公式對能量密度的估計會出現(xiàn)較大的誤差，當(dāng)阻尼很小時，能量密度分布趨向于平均，但近場項等其它因素對結(jié)果的影響會更加明顯。仿真結(jié)果也表明，對于點載荷作用的情況，通過分別計算直接場和混響場能量的方法可以改善EFEA 的計算結(jié)果。