耦合板結(jié)構(gòu)的非結(jié)構(gòu)零階能量有限元分析

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

(中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系,中國(guó)科學(xué)院 材料力學(xué)行為和設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，合肥　230026)

?

耦合板結(jié)構(gòu)的非結(jié)構(gòu)零階能量有限元分析

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

(中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系,中國(guó)科學(xué)院 材料力學(xué)行為和設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，合肥230026)

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

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

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

1非結(jié)構(gòu)零階能量有限元法

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

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

(1)

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

1.2非結(jié)構(gòu)網(wǎng)格有限體積法離散

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

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

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

(2)

由Gauss散度定理，有：

(3)

(4)

(5)

(6)

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

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

∫ViηωedV=ηωeVi

(7)

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

∫Viπin,idV=Πin,i

(8)

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

(9)

其中:

(10)

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

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

Bitsie[12]推導(dǎo)出在無(wú)能量損耗的耦合邊界上，能量流與各板在耦合處的能量密度之間的關(guān)系如下

(11)

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

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

(12)

(13)

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

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

首先將式(9)表示的單元?jiǎng)偠确匠谭謩e代入三種行波的控制方程式(1)，可以得出不含耦合信息的系統(tǒng)總體方程：

(14)

式中:

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

(15)

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

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

(16)

2算例

2.1L型耦合板

L型耦合板(如圖3所示)的材料均為鋁：密度為2 700 kg/m3，彈性模量為71 GPa，泊松比為0.3。假設(shè)三種波場(chǎng)全分析頻率范圍內(nèi)的內(nèi)損耗因子均為0.05。兩板尺寸均為1 m×1 m，板1厚度為5 mm，板2厚度為1 mm。在板1中心處加載1 W的彎曲波能量，加載頻率范圍為1 kHz～10 kHz的11個(gè)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倍頻程N(yùn)o.9(fc=6.3 kHz)內(nèi)各波場(chǎng)的能量密度分布。各波場(chǎng)能量密度在兩板的耦合邊界處均不連續(xù)。對(duì)于彎曲波能量密度，加載點(diǎn)處最高，距加載點(diǎn)越遠(yuǎn)則越?。欢鎯?nèi)波能量密度，耦合邊界處最高，板2中面內(nèi)波能量比板1大。這是由于該結(jié)構(gòu)的外界能量輸入只有板1中心的彎曲波能量，而面內(nèi)波是來(lái)自于彎曲波傳播到耦合邊界上時(shí)的反射和折射,且面內(nèi)波對(duì)彎曲波的透射系數(shù)遠(yuǎn)大于反射系數(shù)。

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

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

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

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

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

通過(guò)對(duì)比，能說(shuō)明本文處理的正確性。另外，相比于SEA，能量有限元能得到如圖5～6所示的空間分布的能量結(jié)果，有助于工程上找到結(jié)構(gòu)的危險(xiǎn)位置。

2.2簡(jiǎn)化的汽車外殼

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

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

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

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

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

圖12　底板的各波場(chǎng)總能量uEFEM0及SEA結(jié)果對(duì)比Fig.12 Comparison of total energy result of bottom plate, by uEFEM0 and SEA

3結(jié)論

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

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

參 考 文 獻(xiàn)

[1]伍先俊, 朱石堅(jiān), 曹建華. 結(jié)構(gòu)聲振研究的功率流方法[J]. 力學(xué)進(jìn)展, 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]孫麗萍, 聶武. 桿, 梁結(jié)構(gòu)振動(dòng)能量密度的簡(jiǎn)化解法[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 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ìn), 李鴻光, 孟光. 耦合板結(jié)構(gòu)隨機(jī)能量有限元分析[J]. 振動(dòng)與沖擊, 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] 劉學(xué)哲, 余云龍, 王瑞利, 等. 非結(jié)構(gòu)任意多邊形網(wǎng)格輻射擴(kuò)散方程有限體積格式[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用, 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月生

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

關(guān)鍵詞：非結(jié)構(gòu)零階能量有限元；耦合板結(jié)構(gòu)；多波場(chǎng)高頻振動(dòng)

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

中圖分類號(hào):TU311.3

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.13.024

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

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