能量有限元方法研究進展及其非平面波情況下的改進

于開平 王懷志

(哈爾濱工業大學航天學院，哈爾濱 150001)

0 引言

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

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

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

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

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

圖1 整流罩結構EFEA 模型及仿真結果Fig.1 EFEA model and simulation results of the fairing

圖2 聲腔的聲壓級響應Fig.2 Sound pressure level response of the acoustic cavity

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

1 能量密度的控制方程

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

圖3 單元體的能量流關系Fig.3 Energy flow relationship of unit element

這個關系可以用以下公式來表示[5]

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

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

其中，η是內損耗因子，ω是圓頻率。

結構中單個面波分量可以用以下公式來描述

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

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

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

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

其中，Cg是彎曲波對應的群速度。

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

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

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

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

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

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

公式(8)可通過Galerkin 加權余量方法來構建能量有限元方程，對于2 維板結構，余量方程為

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

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

其中，

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

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

其中，

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

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

對于公式中的能量輸入項Pin，一般可以通過無限大結構的原點輸入導納來得到，即

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

2 EFEA 結果的有效性

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

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

Gur[33]通過對板梁結構在特定頻帶內響應的研究提出了板和梁結構的波長準則，對于梁結構，一般要求單個梁的長度至少大于5 倍的彎曲波長；對于板結構，板的特征長度要求大于2.4 倍的彎曲波長。

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

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

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

3 非平面波假設的改進

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

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

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

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

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

對于單板系統，直接場能量在到達結構邊界的時候會發生反射，而反射的功率流則作為混響場的能量輸入，即在結構的邊界，能量流滿足公式(22)

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

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

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

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

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

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

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

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

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

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

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

滯后阻尼系數為0.1 時，直接場的能量占主導，此時在遠離激勵點位置的能量以及混響場的能量幾乎可以忽略。基于平面波的公式在這種情況下的結果相對SEA 具有一定的優勢。對于耦合板，不同阻尼的情況下，通過單獨計算直接場能量和混響場能量的方法結果與有限元結果更加接近，這與四邊耦合單板的結果時一致的。

4 結論

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