周紅衛(wèi)， 陳海波， 王用巖

(中國科學技術大學 近代力學系,中國科學院 材料力學行為和設計重點實驗室，合肥　230026)

?

耦合板結構的非結構零階能量有限元分析

周紅衛(wèi)， 陳海波， 王用巖

(中國科學技術大學 近代力學系,中國科學院 材料力學行為和設計重點實驗室，合肥230026)

在模擬結構高頻振動時，隨著分析頻率的增高，傳統(tǒng)的數(shù)值分析方法，如有限元(FEA)和邊界元(BEM)等，面臨著計算成本急劇增大的瓶頸[1]。于是，以統(tǒng)計能量分析(Statistical Energy Analysis，SEA)和能量流分析法(Energy Flow Analysis，EFA)等為代表的高頻分析方法已經(jīng)越來越廣泛地應用于結構的高頻振動及噪聲分析中。

目前廣泛使用的SEA基于模態(tài)理論，主要研究變量是耦合結構子系統(tǒng)的總能量；而EFA基于波動理論，其研究的主變量為波長和周期上雙重平均的能量密度，在中高頻內(nèi)能較為準確地反映能量響應水平在空間上的分布，相比于SEA能得到更為豐富的結果，具有一定的優(yōu)越性。Nefske等[2]利用有限元技術求解了梁內(nèi)彎曲波的EFA偏微分方程，這種求解格式被稱為能量有限元(Energy Finite Element Method，EFEM)。后續(xù)不少學者將EFEM擴展到桿、梁[3-4]、板[5-6]等常見基本結構中。2000年Wang[7]提出了零階能量有限元法(EFEM0)，這種方法結合了SEA和EFEM的優(yōu)點[8]，它利用有限體積法離散EFA偏微分方程，最終可以得到與SEA類似的線性方程組。但是，迄今EFEM0均基于結構網(wǎng)格，需要在板結構上劃分結構化的四邊形網(wǎng)格，這在外形不規(guī)則時往往會出現(xiàn)畸形單元，影響計算精度。如以三角形劃分幾何外形，即可克服上述問題。因此，發(fā)展非結構網(wǎng)格的零階能量有限元法(Unstructured Zero-Order EFEM，uEFEM0)顯得非常有實際意義。

基于以上原因，本文推導了三角形網(wǎng)格的uEFEM0的離散格式。對于耦合結構，結合考慮多波場時計算耦合板結構連接處耦合矩陣的一般處理方法[9]，對L型耦合板及簡化汽車外殼的彎曲、拉伸和剪切波場的能量響應進行了數(shù)值模擬，并利用EFEM及商業(yè)軟件AutoSEA對相同結構分別進行了分析，與本文計算結果進行了對比驗證。

1非結構零階能量有限元法

1.1薄板中的能量流分析方程

在高頻振動的薄板中，存在三種不同形式的行波：彎曲波、拉伸波和剪切波。不同頻率和不同種類的行波都互不相關地傳播，并且它們的能量密度分布滿足以下方程[10]：

(1)

式中:e是空間和時間上雙重平均的能量密度，πin是單位面積上的輸入能量，η是阻尼損耗因子，cg是板中行波的群速度，ω為圓頻率。下標B、L和S分別代表彎曲波、拉伸波和剪切波。

1.2非結構網(wǎng)格有限體積法離散

劉學哲等[11]推導了輻射擴散方程的非結構有限體積法格式。本文基于此，推導了式(1)的離散格式。式(1)中的三種行波的控制方程形式一致，下列推導過程中略去了下標B、L或S。

圖1　相鄰控制體Fig.1 Adjacent control volume

將控制方程(1)在控制體i上進行積分：

(2)

由Gauss散度定理，有：

(3)

(4)

(5)

(6)

圖2　共同邊界上定義的向量Fig.2 Vectors on the common boundary

另外，認為e在Vi內(nèi)恒定，則有：

∫ViηωedV=ηωeVi

(7)

并在控制體i內(nèi)積分輸入能量為：

∫Viπin,idV=Πin,i

(8)

綜合式(5)、(7)和(8)，可得uEFEM0在控制體i上的單元方程：

(9)

其中:

(10)

可以驗證，如果劃分規(guī)則的四邊形網(wǎng)格，由uEFEM0得出的計算格式(9)和Wang[7]推導出的應用結構化網(wǎng)格的EFEM0的計算格式是完全一致的。

1.3耦合邊界上的能量流分析

Bitsie[12]推導出在無能量損耗的耦合邊界上，能量流與各板在耦合處的能量密度之間的關系如下

(11)

式中:qcp為耦合邊界上的能量流向量，I為單位矩陣，τ為功率傳輸系數(shù)矩陣，cg為由群速度組成的對角陣，ecp為各板在耦合邊界上的能量密度。

若考慮編號為i和j的兩塊板耦合，并考慮三種行波的情形，則有如下具體形式：

(12)

(13)

由于不同子系統(tǒng)耦合節(jié)點上的能量密度是不相等的，因而在有限元離散的EFEM模型中需要增加耦合邊界上的額外節(jié)點，稱為重節(jié)點。而在EFEM0中，能量密度被設置在單元中心處，耦合邊界上并無節(jié)點，因此無須設置重節(jié)點，而只需將qj在耦合邊界上的計算格式相應地從式(5)改為式(11)，而后組裝到總體剛度方程即可。由于這一原因，只需找出FEM模型中的耦合邊界，便可將此模型直接用于EFEM0的后續(xù)分析。

1.4總體剛度方程的具體形式

首先將式(9)表示的單元剛度方程分別代入三種行波的控制方程式(1)，可以得出不含耦合信息的系統(tǒng)總體方程：

(14)

式中:

eP={e1Pe2P…eNP}，P=B,L式S

(15)

N為系統(tǒng)的控制體個數(shù)，BB等為相應波形的總體剛度矩陣。

接著找出系統(tǒng)中所有的耦合單元對，并應用式(11)添加耦合邊界處各波場能量相互轉(zhuǎn)化對于總體剛度矩陣的影響，得到最終的總體剛度矩陣為：

(16)

2算例

2.1L型耦合板

L型耦合板(如圖3所示)的材料均為鋁：密度為2 700 kg/m3，彈性模量為71 GPa，泊松比為0.3。假設三種波場全分析頻率范圍內(nèi)的內(nèi)損耗因子均為0.05。兩板尺寸均為1 m×1 m，板1厚度為5 mm，板2厚度為1 mm。在板1中心處加載1 W的彎曲波能量，加載頻率范圍為1 kHz～10 kHz的11個1/3倍頻程，各1/3倍頻程中心頻率fc如表1所示。用商業(yè)軟件ANSYS劃分了如圖4所示的三套網(wǎng)格，稱為Mesh1～Mesh3，各模型的單元數(shù)分別為100、400和444。

表1　各1/3倍頻程的中心頻率

圖3　彎曲波輸入的L型耦合板Fig.3 L-shape plates with bending wave load

圖4　L型板的三套網(wǎng)格Fig.4 Three mesh models of L-shape plates

圖5顯示了Mesh3在1/3倍頻程No.9(fc=6.3 kHz)內(nèi)各波場的能量密度分布。各波場能量密度在兩板的耦合邊界處均不連續(xù)。對于彎曲波能量密度，加載點處最高，距加載點越遠則越??；而面內(nèi)波能量密度，耦合邊界處最高，板2中面內(nèi)波能量比板1大。這是由于該結構的外界能量輸入只有板1中心的彎曲波能量，而面內(nèi)波是來自于彎曲波傳播到耦合邊界上時的反射和折射,且面內(nèi)波對彎曲波的透射系數(shù)遠大于反射系數(shù)。

圖5　耦合板各波場的能量響應；參考值：1x10-12 J/m2Fig.5 Energy density of L-shape plates; ref=1x10-12 J/m2

圖6中為三種網(wǎng)格下的相同位置(如圖3中的虛線Line1,1/3倍頻程No.11)各波場能量密度的對比情況，橫坐標為Line1上距點A的距離，點A、點B和點C為各邊的中點?？梢钥闯觯瑄EFEM0具有良好的網(wǎng)格收斂性。另外，在亂序網(wǎng)格Mesh3上良好的適應性也說明了uEFEM0能適用于不規(guī)則結構的三角形網(wǎng)格劃分。Park等[13]通過算例指出，隨著分析頻率的增高和內(nèi)損耗因子的增大，在板2末端的面內(nèi)波能量密度可能超過彎曲波能量密度，圖6取得了與之類似的結果。這就使得在阻尼偏大或頻率較高的情況下，考察距離加載點較遠的振動結構能量響應時，面內(nèi)波往往不能忽略。

圖6　三套網(wǎng)格的相同位置的能量密度結果對比Fig.6 Comparison of energy density result of three mesh models on the same position

另外，對該結構用常規(guī)EFEM(劃分如Mesh2網(wǎng)格相同尺寸的四邊形網(wǎng)格，僅計算了彎曲波)以及商業(yè)軟件AutoSEA分別進行了分析。三種方法得到的各板總能量結果對比情況(Mesh2)如圖7～9所示。從各圖中可以看出，EFEM和uEFEM0的彎曲波計算結果基本完全吻合，這說明本文提出的uEFEM0計算精度與傳統(tǒng)EFEM一致。另外，uEFEM0和AutoSEA得到的各波場結果吻合較好。板2的彎曲波能量有1.5 dB左右的差距，板1的高頻拉伸波結果有1～1.5dB左右的差距。這是由于能量有限元和統(tǒng)計能量法兩種方法的差異造成的。相比較游進等得到的結果[6]，本文考慮了全部三種波形，拉伸波的對比結果更好一些。

圖7 L型耦合板彎曲波總能量uEFEM0、EFEM和SEA結果對比Fig.7ComparisonoftotalbendingenergyresultofL-shapeplates,byuEFEM0,EFEMandSEA圖8 L型耦合板拉伸波總能量uEFEM0和SEA結果對比Fig.8ComparisonoftotallongitudinalenergyresultofL-shapeplates,byuEFEM0andSEA圖9 L型耦合板剪切波總能量uEFEM0和SEA結果對比Fig.9ComparisonoftotalshearenergyresultofL-shapeplates,byuEFEM0andSEA

通過對比，能說明本文處理的正確性。另外，相比于SEA，能量有限元能得到如圖5～6所示的空間分布的能量結果，有助于工程上找到結構的危險位置。

2.2簡化的汽車外殼

圖10　簡化的汽車外殼的外形及其uEFEM0模型Fig.10 uEFEM0 model of simplified vehicle shell

如圖10所示，簡化的汽車外殼材料為鋁(參數(shù)同L型耦合板算例)，各板厚度均為5 mm。uEFEM0模型的網(wǎng)格數(shù)為468，結構中共有耦合單元79對。分析的頻率范圍為1 kHz～10 kHz(如表1)，假設全頻域內(nèi)的各波場阻尼損耗因子均為0.06。車體尺寸為：車身長度3 m，寬度1.5 m，高度1 m；車頭高度0.4 m，頂蓋長度2 m。設4個車輪在底板坐標系XOY內(nèi)的坐標為(0.6,0.3)、(0.6,1.2)、(2.5,0.3)和(2.5,1.2)。該汽車殼側板形狀不規(guī)則，用結構網(wǎng)格不易劃分,但易于劃分非結構化網(wǎng)格。該案例可以初步體現(xiàn)非結構網(wǎng)格在處理不規(guī)則形狀時的優(yōu)勢。

在車身底板4個車輪處施加彎曲波能量，假設分析頻域內(nèi)的4個外部功率輸入均為0.25 W。圖11從兩個角度顯示了結構在中心頻率為fc=6.3 kHz的1/3倍頻程內(nèi)的彎曲波能量密度分布。圖中的單位為dB，參考值為1x10-12J/m2?？梢钥闯?，底板上加載激勵的車輪處的能量最高，而據(jù)底板最遠的頂板處能量最低。

圖11　簡易汽車外殼的彎曲波能量密度分布；fc =6.3 kHzFig.11 Bending energy density distribution of simplified vehicle shell. fc =6.3 kHz

用AutoSEA建立對應的SEA模型，并在底板子系統(tǒng)全分析頻段內(nèi)加載彎曲波能量1 W。圖12為全分析頻段內(nèi)汽車底板的總能量結果對比。從對比結果中可以看出，兩種方法的彎曲波能量結果很接近，面內(nèi)波能量有大約1～2 dB的差距，這和L型耦合板算例內(nèi)取得的結果類似。

圖12　底板的各波場總能量uEFEM0及SEA結果對比Fig.12 Comparison of total energy result of bottom plate, by uEFEM0 and SEA

3結論

本文基于薄板結構中三種波場EFA控制方程，推導了應用三角形網(wǎng)格的非結構網(wǎng)格零階能量有限元的離散格式，提出了用于求解薄板高頻振動的非結構零階能量有限元(uEFEM0)分析方法。并運用uEFEM0計算了兩個算例，均取得了與EFEM和AutoSEA較為一致的結果，驗證了本文提出的方法的正確性。零階能量有限元在耦合邊界上不需要設置重節(jié)點，簡化了使用FEM模型進行高頻能量流分析的步驟。相比于EFEM0，由于uEFEM0采用三角形網(wǎng)格，能適應更為復雜的結構外形，具有更為廣闊的工程應用前景。

另外考慮到三維聲空間易于用四面體網(wǎng)格劃分，使得本文方法在考慮聲振耦合方面具有較大的優(yōu)勢。

參 考 文 獻

[1]伍先俊, 朱石堅, 曹建華. 結構聲振研究的功率流方法[J]. 力學進展, 2006, 36(3): 363-372.

WU Xian-jun, ZHU Shi-jian, CAO Jian-hua. Review on structural vibration power flow prediction methods[J]. Advances in Mechanics, 2006, 36(3): 363-372.

[2]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 Stress and Reliability in Design, 1989, 111(1): 94-100.

[3]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.

[4]孫麗萍, 聶武. 桿, 梁結構振動能量密度的簡化解法[J]. 哈爾濱工程大學學報, 2004, 25(4): 403-406.

SUN Li-ping, NIE Wu. Asimplified method for solving energy density of rods and beams[J]. Journal of Harbin Engineering University, 2004, 25(4):403-406.

[5]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-164.

[6]游進, 李鴻光, 孟光. 耦合板結構隨機能量有限元分析[J]. 振動與沖擊, 2009, 28(11): 43-46.

YOU Jin, LI Hong-guang, MENG Guang. Random energy finite element analysis of coupled plate structures[J]. Journal of Vibration and Shock, 2009, 28(11):43-46.

[7]Wang Shuo. High frequency energy flow analysis methods: Numerical implementation, applications, and verification[D]. Purdue University: United States, 2000.

[8]Moens I. On the use and the validity of the energy finite element method for high frequency vibrations[D]. Catholic University of Louvain,2001.

[9]Langley R S, Heron K H. Elastic wave transmission through plate/beam junctions[J]. Journal of Sound and Vibration, 1990, 143(2): 241-253.

[10] Dong J, Choi K K, Wang A, et al. Parametric design sensitivity analysis of high‐frequency structural-acoustic problems using energy finite element method[J]. International Journal for Numerical Methods in Engineering, 2005, 62(1): 83-121.

[11] 劉學哲, 余云龍, 王瑞利, 等. 非結構任意多邊形網(wǎng)格輻射擴散方程有限體積格式[J]. 數(shù)值計算與計算機應用, 2010,31(4): 259-270.

LIU Xue-zhe, YU Yun-long, WANG Rui-li, et al. A cell-centered finite volume scheme for discretizing diffusion equation on unstructured arbitrary polygonal meshes[J], Journal on Numerical Methods and Computer Applications, 2010, 31(4): 259-270.

[12] Bitsie F. The structural-acoustic energy finite-element method and energy boundary-element method[D]. Purdue University: United States, 1996.

第一作者 周紅衛(wèi) 男，博士生，1987年6月生

摘要：基于薄板高頻振動時三種波場的能量流控制方程，用三角形網(wǎng)格劃分板結構，推導了非結構零階能量有限元的計算格式。并運用該方法計算了L型耦合板和簡化的汽車外殼，得到各波場在空間上的能量密度結果。將計算結果與統(tǒng)計能量法和能量有限元進行對比,驗證了其正確性。計算結果表明，頻率很高時面內(nèi)波能量往往能達到彎曲波能量水平，高頻振動仿真時有必要同時考慮三種波場。

關鍵詞：非結構零階能量有限元；耦合板結構；多波場高頻振動

Unstructured zero-order energy finite element method for coupled plate structures

ZHOUHong-wei,CHENHai-bo,WANGYong-yan(CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, Hefei 230026, China)

Abstract:Based on the governing equations of energy flow analysis (EFA), meshing plates with triangles, an unstructured zero-order energy finite element method (uEFEM0) was developed. A procedure used to calculate L-shape coupled plates and a simplified vehicle shell was presented considering both bending and in-plane wave fields. The proposed method was used to predict the distribution of energy response of plates. To confirm its validity, the energy finite element method (EFEM) and statistical energy analysis (SEA) were employed to simulate the same structures, and the results agreed well with those using the proposed method. It was shown that in order to simulate plates vibration within a higher frequency range, it is necessary to consider not only the bending wave field but also the in-plane wave field, since the in-plane wave energy level may be close to the bending wave energy level.

Key words:unstructured zero-order finite element; coupled plates structure; high frequency vibration of multi-wave fields

中圖分類號:TU311.3

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.13.024

通信作者陳海波 男，教授，1968年2月生

收稿日期：2014-03-03修改稿收到日期：2014-07-23